10.1 Introdução


Em um experimento de superfície de resposta, as variáveis ou fatores independentes podem variar em um intervalo contínuo. O objetivo é determinar as configurações dos fatores que produzem uma resposta máxima ou mínima ou mapear a relação entre a resposta e as configurações dos fatores neste espaço de fatores contíguos.

Do ponto de vista prático, se quisermos saber muito sobre a relação entre os factores e a resposta, serão necessárias muitas experiências. Por essa razão, os projetos de superfícies de resposta raramente são conduzidos com mais de seis fatores. Experimentos de superfície de resposta são normalmente usados no último estágio de experimentação. Os factores importantes já foram determinados em experiências anteriores e nesta fase da experimentação o objectivo é descrever em detalhe a relação entre os factores e a resposta.

Geralmente é conhecido ou assumido que um modelo linear simples, mesmo com interações, não é bom o suficiente para representar essa relação. Para localizar máximos ou mínimos na resposta em função das configurações do fator, pelo menos três níveis de cada fator devem ser utilizados.

Os métodos de superfície de resposta geralmente se referem a um pacote completo de ferramentas estatísticas de planejamento e análise que são usadas nas três etapas a seguir.

Os métodos de superfície de resposta têm encontrado uso considerável na indústria, especialmente em processos químicos onde o rendimento da reação ou o custo de produção podem ser otimizados em função das configurações de fatores de processo controláveis. Desde a sua origem, esses planos também encontraram aplicações bem-sucedidas em ciência de alimentos, engenharia, biologia, psicologia, têxteis, educação e muitas outras áreas.


10.2 Fundamentos da metodologia de superfície de resposta



10.2.1 Modelo quadrático empírico


Nos métodos de superfície de resposta, a relação entre a resposta \(y\) e as configurações do fator \(x\)’s é assumida como uma equação não linear dada por \[ y = f(x)+\epsilon\cdot \]

Para dois fatores ou variáveis independentes, poderíamos escrever esta equação como \[\begin{equation*} \tag{10.1} y = f (x_1, x_2) + \epsilon, \end{equation*}\] onde \(f\) é uma função não linear e \(\epsilon\) é o efeito aleatório do erro experimental. Quando \(f\) é desconhecido, pode ser aproximado perto do ponto \((x_{10}, x_{20})\) usando a aproximação em Série de Taylor bi-dimensional, ou seja, \[\begin{equation*} \tag{10.2} \begin{array}{rcl} f(x_1,x_2) & \approx & f(x_{10},x_{20})+(x_1-x_{10})\left. \dfrac{\partial f(x_1,x_2)}{\partial x_1}\right|_{x_1=x_{10},x_2=x_{20}} \\ & & \qquad \qquad +(x_2-x_{20})\left. \dfrac{\partial f(x_1,x_2)}{\partial x_1}\right|_{x_1=x_{10},x_2=x_{20}} \\ & & \qquad \qquad +\dfrac{1}{2}(x_1-x_{10})^2\left. \dfrac{\partial^2 f(x_1,x_2)}{\partial x_1^2}\right|_{x_1=x_{10},x_2=x_{20}} \\ & & \qquad \qquad +\dfrac{1}{2}(x_2-x_{20})^2\left. \dfrac{\partial^2 f(x_1,x_2)}{\partial x_2^2}\right|_{x_1=x_{10},x_2=x_{20}} \\ & & \qquad \qquad +\dfrac{1}{2}(x_1-x_{10})(x_2-x_{20})\left. \dfrac{\partial^2 f(x_1,x_2)}{\partial x_1\partial x_2}\right|_{x_1=x_{10},x_2=x_{20}}, \\ \end{array} \end{equation*}\] o que leva a uma equação quadrática geral da forma \[\begin{equation*} \tag{10.3} y=\beta_0+\beta_1 x_1+\beta_2 x_2+\beta_{11}x_1^2 +\beta_{22}x_2^2+\beta_{12}x_1x_2+\epsilon \end{equation*}\] onde \[ \beta_1=\left. \dfrac{\partial f(x_1,x_2)}{\partial x_1}\right|_{x_1=x_{10},x_2=x_{20}} \] e assim por diante.

Se a região de interesse for de tamanho moderado, esta equação quadrática geral fornecerá uma boa aproximação para \(f\) e pode ser usada para interpolação dentro desta região.

A equação quadrática geral é bastante flexível e com coeficientes apropriados ela pode descrever uma ampla variedade de superfícies, como topos de colinas, vales, cristas ascendentes ou descendentes, ou pontos de sela, como mostrado nos gráficos de contorno na Figura 10.1.

Figura 10.1: Superfícies que podem ser descritas pela equação quadrática geral.

Com \(k\) fatores ou variáveis independentes, a equação quadrática geral pode ser escrita na forma \[\begin{equation*} \tag{10.4} \end{equation*}\] e a menos que a função \(f\) seja conhecida, esta equação forma a base dos métodos de superfície de resposta.

Projetos de superfície de resposta foram criados para fornecer dados para aproximar esta equação e ferramentas matemáticas foram criadas para explorar a superfície ajustada representada por esta equação.


10.2.2 Considerações de planejamento


O primeiro requisito de um planejamento de superfície de resposta é que ele forneça dados que permitirão a estimação dos coeficientes no modelo (10.4). Este modelo possui \(1+2k+k(k-1)/2\) coeficientes, portanto, qualquer projeto de superfície de resposta deve ter pelo menos três níveis para cada fator, isto para permitir a estimação de termos quadráticos e pelo menos \(1+2k +k(k−1)/2\) execuções totais.

Os experimentos fatoriais \(3^k\) têm três níveis para cada fator. No entanto, Box and Wilson (1951) mostraram que eles eram menos satisfatórios como um projeto de superfície de resposta do que um projeto alternativo que eles chamaram de projeto composto central. O projeto composto central será o primeiro projeto de superfície de resposta padrão que discutiremos na próxima seção.

Normalmente os fatores ou variáveis independentes \(x_i\) na Equação quadrática geral (10.4) são codificados e escalonados conforme mostrado nas Seções 3.7.2 e 3.7.3 de modo que a região experimental seja um hipercubo ou hiperesfera com raio \(R\) no espaço codificado.

Este modelo quadrático geral pode ser escrito em termos matriciais como \[\begin{equation*} \tag{10.5} y=x\beta+x^\top\pmb{B}x+\epsilon, \end{equation*}\] onde \(x^\top=(1,x_1,x_2,\cdots,x_k)\), \(\beta^\top=(\beta_0,\beta_1,\cdots,\beta_k)\) e \(\pmb{B}\) é a matriz simétrica \[ \pmb{B}=\begin{pmatrix} \beta_{11} & \beta_{12}/2 & \cdots & \beta_{1k}/2 \\ & \beta_{22} & \cdots & \beta_{2k}/2 \\ & & \ddots & \\ & & & \beta_{kk} \end{pmatrix} \]

Ao ajustar um modelo de regressão linear da forma \(y = x\beta\), os pontos do projeto são escolhidos para minimizar a variância dos coeficientes ajustados \(\widehat{\beta} = (X^\top X)^{-1} Xy\). Como a matriz de variância e covariância dos coeficientes de regressão estimados é \(\sigma^2(X^\top X)^{-1}\), isso significa que os pontos do projeto devem ser escolhidos de forma que a matriz \((X^\top X)\) seja diagonal ou que o projeto seja ortogonal como os projetos ou planejamentos \(2^k\) discutido no Capítulo 3 e os projetos \(2^{k-p}\) discutidos no Capítulo 6.

Quando a matriz \((X^\top X)\) é diagonal, os elementos diagonais de \((X^\top X)^{-1}\) serão minimizados. Por outro lado, ao ajustar o modelo quadrático geral, o objetivo principal não é compreender o mecanismo da relação subjacente entre a resposta e os fatores.

Portanto, os coeficientes específicos no modelo quadrático geral são de menor importância. O que é mais importante em um estudo de superfície de resposta é desenvolver uma equação de previsão com o objetivo final de determinar as condições operacionais ideais. Assim, a variância de um valor previsto em \(x\), que é dada pela equação \[ \mbox{Var}\big(\widehat{y}(x)\big) = \sigma^2 x^\top (X^\top X)^{-1}x \] é mais importante do que a variância dos coeficientes ajustados.

Como não se sabe antes do experimento ser conduzido onde o ótimo estará na região de projeto, uma propriedade desejável de um projeto de superfície de resposta é ter a variância de um valor previsto quase igual em todos os lugares da região de projeto. Desta forma, a precisão na previsão do ótimo não dependerá da sua posição desconhecida na região de projeto.

O primeiro passo para equalizar a variância da previsão na região do projeto é encontrar um projeto que seja rotacionável (Myers and Montgomery, 2002). Um projeto rotacionável é aquele em que a variância do valor previsto no ponto \(x\) é apenas uma função da distância da origem do projeto até \(x\). O projeto composto central de Box and Wilson e outros projetos de superfície de resposta padrão têm essa propriedade.

Box and Hunter (1957) mostraram que certos projetos rotacionáveis poderiam ser modificados para ter precisão uniforme, o que significa que a variância de um valor previsto é a mesma na origem e no raio de um na região de projeto codificada.

Um projeto com precisão uniforme está próximo de ter a variância de um valor previsto igual em toda a região de projeto. Muitos projetos de superfície de resposta padrão possuem essa propriedade ou estão próximos de possuí-la. Quando um projeto de superfície de resposta padrão não atende às necessidades de um experimentador, um algoritmo de computador pode ser usado para construir um projeto de propósito especial.

Seria razoável usar um projeto que fosse D-ótimo, descrito pela primeira vez na Seção 6.5.2, para o modelo (10.4) como um projeto de superfície de resposta. No entanto, o critério D-otimalidade minimiza a variância dos coeficientes do modelo. O critério de I-otimalidade minimiza a variância média de um valor previsto e, portanto, é mais apropriado para encontrar um projeto de superfície de resposta.

Programas, como a função optFederov no pacote R AlgDesign, podem ser usados para criar planejamentos I-ótimos e D-ótimos. Giovannitti-Jensen and Myers (1989) e Myers et al. (1992) descrevem um gráfico de dispersão de variância que pode ser usado para avaliar visualmente se qualquer projeto de superfície de resposta está próximo de ter as propriedades de rotacionabilidade e precisão uniforme.

Além de fornecer uma boa distribuição da variância de um valor previsto, duas outras propriedades desejáveis dos projetos de superfície de resposta são

  1. fornecer uma estimativa do erro experimental “puro” para que a adequação do modelo quadrático geral possa ser verificada, e

  2. permitir o bloqueio para que um experimentador possa começar com um projeto linear e adicionar um segundo bloco para estimar a curvatura, se necessário.

Para estimar o erro puro, pelo menos um ponto do projeto deve ser replicado. Ao bloquear, o primeiro bloco seria consistem em um projeto como \(2^k\) ou \(2^{k−p}\) aumentado com pontos centrais que permitem verificar a adequação do modelo linear.

Se o modelo linear for adequado, nenhuma experimentação adicional será necessária. Caso o modelo linear não seja adequado, um segundo bloco de experimentos poderá ser adicionado que permitirá estimar os coeficientes quadráticos no modelo quadrático geral.


10.3 Projetos padrão para modelos de segunda ordem


Esta seção apresenta alguns dos designs de superfície de resposta completamente aleatória (CRRS) mais populares.


10.3.1 Projeto composto central (CCD)


O planejamento composto central (CCD) de Box and Wilson (1951) consiste em um planejamento fatorial \(2^k\) ou fatorial fracionário \(2^{k-p}\) aumentado com pontos centrais e pontos axiais, como mostrado na Figura 10.2. O design \(2^k\) ou resolução V \(2^{k-p}\) permite a estimativa dos termos lineares e lineares por lineares no modelo quadrático geral.

Figura 10.2: Planejamento composto central em duas e três dimensões.

A adição dos pontos centrais e pontos axiais permite estimar dos termos quadráticos. Ao escolher a distância da origem aos pontos axiais, \(\alpha\) em unidades codificadas, igual a \(\sqrt[4]{F}\) onde \(F\) é o número de pontos na porção fatorial do dimensionamento, um projeto composto central será rotacionável. Ao escolher o número correto de pontos centrais, o projeto composto central terá a propriedade de precisão uniforme.

Como exemplo de um projeto composto central (CCD), considere os resultados do experimento descrito por Kamon et al. (1999) na Tabela 10.1. O experimento foi conduzido para encontrar a formulação ideal de um aditivo de dois componentes para melhorar a trabalhabilidade de argamassas de cimento.

\[ \begin{array}{cccccccc}\hline \mbox{Corrida} & x_1 & x_2 & x_3 & \mbox{Água/Cimento} & \mbox{Licor negro} & \mbox{SNF} & y \\ \hline 1 & -1 & -1 & -1 & 0.330 & 0.120 & 0.080 & 109.5 \\ 2 & 1 & -1 & -1 & 0.350 & 0.120 & 0.080 & 120.0 \\ 3 & -1 & 1 & -1 & 0.330 & 0.180 & 0.080 & 110.5 \\ 4 & 1 & 1 & -1 & 0.350 & 0.180 & 0.080 & 124.5 \\ 5 & -1 & -1 & 1 & 0.330 & 0.120 & 0.120 & 117.0 \\ 6 & 1 & -1 & 1 & 0.350 & 0.120 & 0.120 & 130.0 \\ 7 & -1 & 1 & 1 & 0.330 & 0.180 & 0.120 & 121.0 \\ 8 & 1 & 1 & 1 & 0.350 & 0.180 & 0.120 & 132.0 \\ 9 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 117.0 \\ 10 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 117.0 \\ 11 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 115.0 \\ 12 & -1.68 & 0 & 0 & 0.323 & 0.150 & 0.100 & 109.5 \\ 13 & 1.68 & 0 & 0 & 0.357 & 0.150 & 0.100 & 134.0 \\ 14 & 0 & -1.68 & 0 & 0.340 & 0.100 & 0.100 & 120.0 \\ 15 & 0 & 1.68 & 0 & 0.340 & 0.200 & 0.100 & 121.0 \\ 16 & 0 & 0 & -1.68 & 0.340 & 0.150 & 0.066 & 115.0 \\ 17 & 0 & 0 & 1.68 & 0.340 & 0.150 & 0.134 & 127.0 \\ 18 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 116.0 \\ 19 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 117.0 \\ 20 & 0 & 0 & 0 & 0.340 & 0.150 & 0.100 & 117.0 \\\hline \end{array} \] Tabela 10.1: Projeto composto central (CCD) em unidades codificadas e reais para experimento de trabalhabilidade do cimento.

Os fatores foram a proporção água/cimento, a porcentagem de licor negro adicionado e a porcentagem de SNF adicionado. O lado esquerdo da tabela mostra os fatores codificados \(x_1-x_3\). Os pontos axiais foram escolhidos em \(\pm 1.68 = \sqrt[4]{8}\) em unidades codificadas para tornar o desenho giratório. Licor negro, também conhecido como licor preto, é um subproduto do processo de tratamento químico da indústria de papel e celulose.

Os níveis codificados são encontrados como \[ (\mbox{nível real} - \mbox{valor central})/(\mbox{meio intervalo})\cdot \] Por exemplo, \(x_2 = (\mbox{Licor negro} − 0.150)/0.03\). Os níveis reais dos fatores no lado direito da tabela podem ser obtidos a partir dos valores codificados no lado esquerdo da tabela, resolvendo \[ \mbox{nível real} = (\mbox{meio intervalo})\times x_i + \mbox{valor central}\cdot \]

Ao incluir seis pontos centrais, o plano tem precisão uniforme e a variância de um valor previsto será a mesma na origem \((0, 0, 0)\), em níveis de fator codificados e em qualquer ponto no raio 1 em unidades codificadas. Box and Hunter (1957) tabularam o número de pontos centrais necessários para tornar um projeto composto central com precisão uniforme para vários valores de \(k = \mbox{número de fatores}\).

O código R abaixo chama a função ccd do pacote R rsm, Lenth (2013) e Lenth (2009), para criar um projeto composto central de precisão uniforme em três fatores. As colunas de fator do objeto de planejamento resultante rotd são copiadas no quadro de dados rotdm e a função Vdgraph do pacote R Vdgraph, Lawson (2013b) e Lawson (2012), é então chamada para criar o gráfico de dispersão de variância porposto por Myers et al.’s (1992), para este projeto mostrado na Figura 10.3.

library(rsm)
rotd = ccd(3, n0 = c(4,2), alpha = "rotatable", randomize = FALSE)
rotdm = rotd[ , 3:5]
library(Vdgraph)
par( mar = c(3,3,0,0))
Vdgraph(rotdm)
## number of design points= 20 
## number of factors= 3
##           Radius   Maximum   Minimum   Average
##  [1,] 0.00000000  3.326805  3.326805  3.326805
##  [2,] 0.08660254  3.320828  3.320828  3.320828
##  [3,] 0.17320508  3.303837  3.303837  3.303837
##  [4,] 0.25980762  3.278640  3.278640  3.278640
##  [5,] 0.34641016  3.249923  3.249923  3.249923
##  [6,] 0.43301270  3.224241  3.224241  3.224241
##  [7,] 0.51961524  3.210026  3.210026  3.210026
##  [8,] 0.60621778  3.217583  3.217583  3.217583
##  [9,] 0.69282032  3.259089  3.259089  3.259089
## [10,] 0.77942286  3.348596  3.348596  3.348596
## [11,] 0.86602540  3.502029  3.502029  3.502029
## [12,] 0.95262794  3.737186  3.737186  3.737186
## [13,] 1.03923048  4.073740  4.073740  4.073740
## [14,] 1.12583302  4.533236  4.533236  4.533236
## [15,] 1.21243557  5.139093  5.139093  5.139093
## [16,] 1.29903811  5.916603  5.916603  5.916603
## [17,] 1.38564065  6.892934  6.892934  6.892934
## [18,] 1.47224319  8.097125  8.097125  8.097125
## [19,] 1.55884573  9.560089  9.560089  9.560089
## [20,] 1.64544827 11.314613 11.314613 11.314613
## [21,] 1.73205081 13.395357 13.395357 13.395357
grid()

Figura 10.3: Gráfico de dispersão das variâncias para CCD de precisão uniforme em três fatores.

O gráfico de dispersão de variância traça a variância em escala máxima, mínima e média de um valor previsto \(\mbox{Var}\big(\widehat{y}(x)\big)/\sigma^2\), como uma função da distância da origem do projeto. Quando o desenho não for rotacional, haverá três linhas distintas no gráfico. Quanto mais próximas as linhas, mais próximo o planejamento estará de ser rotacional. Na Figura 10.3 as linhas se sobrepõem mostrando que o desenho é rotacional.

Também pode ser visto que o valor da variância padronizada de um valor previsto na origem é quase igual à variância no raio 1.0, indicando que o projeto tem precisão uniforme.

Os experimentos na Tabela 10.1 devem ser executados em ordem aleatória, mas estão listados em ordem não aleatória para ilustrar outro aspecto útil dos experimentos compostos centrais. Este planejamento pode ser executado como uma sequência de dois blocos. Em algumas aplicações esta é uma abordagem desejável.

As primeiras oito execuções da tabela representam um fatorial \(2^3\) completo. As execuções 9 a 11 são pontos centrais no nível médio de cada fator. Se esses 11 experimentos fossem concluídos no primeiro bloco, os dados poderiam ser analisados para ver se um modelo linear \[\begin{equation*} \tag{10.6} y=\beta_0+\beta_1 x_1 +\beta_2 x_2+\beta_3 x_3+\beta_{12} x_1 x_2 +\beta_{13} x_1 x_3 +\beta_{23} x_2 x_3+\epsilon \end{equation*}\] é adequado para representar os dados.

Nesse caso, as execuções número 12 a 20 não precisam ser concluídas. Se o modelo linear não for adequado para representar os dados, o segundo conjunto de experimentos pode ser concluído e uma variável de bloco adicionada ao modelo para levar em conta quaisquer mudanças nas variáveis de fundo que ocorrem entre o primeiro e o segundo conjunto de experimentos.


10.3.2 Planejamento de Box-Behnken


Enquanto o planejamento composto central requer cinco níveis \((-\alpha, -1, 0, +1, +\alpha)\) para cada fator, Box and Behnken (1960) desenvolveram alguns planejamentos de três níveis que permitirão a estimativa do modelo quadrático geral.

Esses projetos consistem em \(2^2\) fatoriais em cada par de fatores, com todos os outros fatores mantidos constantes em seu nível médio mais alguns pontos centrais. Nenhum projeto Box-Behnken existe por apenas dois fatores. Um exemplo de projeto Box-Behnken em três fatores codificados é mostrado na Tabela 10.2.

\[ \begin{array}{cccccccc}\hline \mbox{Corrida} & x_1 & x_2 & x_3 \\\hline 1 & -1 & -1 & 0 \\ 2 & 1 & -1 & 0 \\ 3 & -1 & 1 & 0 \\ 4 & 1 & 1 & 0 \\ 5 & -1 & 0 & -1 \\ 6 & 1 & 0 & -1 \\ 7 & -1 & 0 & 1 \\ 8 & 1 & 0 & 1 \\ 9 & 0 & -1 & -1 \\ 10 & 0 & 1 & -1 \\ 11 & 0 & -1 & 1 \\ 12 & 0 & 1 & 1 \\ 13 & 0 & 0 & 0 \\ 14 & 0 & 0 & 0 \\ 15 & 0 & 0 & 0 \\\hline \end{array} \] Tabela 10.2: Projeto Box-Behnken em três fatores.

Os projetos Box-Behnken têm duas vantagens sobre os planejamentos compostos centrais. A primeira vantagem é que exigem apenas que os fatores sejam variados em três níveis. Isto pode tornar a experimentação menos dispendiosa se protótipos reais estiverem sendo construídos na experimentação.

A segunda vantagem é que eles geralmente, exceto para o caso de fator com 5 níveis, exigem menos corridas totais do que o planejamento composto central. Por exemplo, o planejamento composto central de três fatores exigiu 20 execuções, enquanto o projeto Box-Behnken de três fatores exigiu apenas 15 execuções.

Uma desvantagem de um projeto Box-Behnken comparado a um planejamento composto central é que ele não pode ser construído em duas etapas, começando com um projeto \(2^k\). Isso pode ser visualizado na Figura 10.4, que diagrama um experimento Box-Behnken de três fatores. Um dimensionamento \(2^3\) consiste em passagens em todos os cantos do cubo, como mostrado na Seção 3.7.1, mas como pode ser visto na Figura 10.4, o dimensionamento Box-Behnken não inclui esses pontos.

Figura 10.4: Representação gráfica do experimento Box-Behnken de três fatores.

Portanto, um projeto Box-Behnken deve ser usado se o experimentador estiver razoavelmente certo de que um modelo linear não representará adequadamente a relação entre os fatores e a resposta, e ele ou ela quiser economizar algumas execuções. Outra possível desvantagem é que os fatores possuem apenas três níveis. Embora isto possa ser menos dispendioso na construção de protótipos, ter três níveis não deixa nada para verificar a adequação do modelo quadrático.

A Figura 10.5 mostra o gráfico de dispersão de variância para o experimento Box-Behnken de três fatores. Como as linhas para a variância escalonada mínima, máxima e média de um valor previsto não se sobrepõem como na Figura 10.3, o projeto Box-Behnken não é rotacionável. Mas, como as linhas estão muito próximas dentro da região experimental codificada com um raio de 1.0, o projeto Box-Behnken é próximo o suficiente para ser rtotacionável e ter precisão uniforme para uso prático. Isto é verdade para projetos Box-Behnken com \(k = 3\) a 6.

# create design with rsm
library(rsm)
bbd3 = bbd(3,randomize=FALSE,n0=3)
library(Vdgraph)
Vdgraph(bbd3[ , 3:5])
## number of design points= 15 
## number of factors= 3
##           Radius   Maximum   Minimum   Average
##  [1,] 0.00000000  5.000000  5.000000  5.000000
##  [2,] 0.08660254  4.984477  4.984445  4.984458
##  [3,] 0.17320508  4.939125  4.938625  4.938825
##  [4,] 0.25980762  4.867602  4.865070  4.866083
##  [5,] 0.34641016  4.776000  4.768000  4.771200
##  [6,] 0.43301270  4.672852  4.653320  4.661133
##  [7,] 0.51961524  4.569125  4.528625  4.544825
##  [8,] 0.60621778  4.478227  4.403195  4.433208
##  [9,] 0.69282032  4.416000  4.288000  4.339200
## [10,] 0.77942286  4.400727  4.195695  4.277708
## [11,] 0.86602540  4.453125  4.140625  4.265625
## [12,] 0.95262794  4.596352  4.138820  4.321833
## [13,] 1.03923048  4.856000  4.208000  4.467200
## [14,] 1.12583302  5.260109  4.367570  4.724583
## [15,] 1.21243557  5.839134  4.638625  5.118825
## [16,] 1.29903811  6.625977  5.043945  5.676758
## [17,] 1.38564065  7.656000  5.608000  6.427200
## [18,] 1.47224319  8.966977  6.356945  7.400958
## [19,] 1.55884573 10.599125  7.318625  8.630825
## [20,] 1.64544827 12.595102  8.522570 10.151583
## [21,] 1.73205081 15.000000 10.000000 12.000000
grid()

Figura 10.5: Gráfico de dispersão de variância para experimento Box-Behnken em três fatores.

Outro gráfico que é útil para avaliar a variância dos valores previstos de um projeto de superfície de resposta é o gráfico da fração do espaço de planejamento (Zahran et al., 2003). Este gráfico mostra a variância relativa de um valor previsto, \(\mbox{Var}\big(\widehat{y}(x)\big)/\sigma^2\), no eixo vertical versus a fração de pontos na região de projeto no eixo horizontal. A Figura 10.6 mostra o gráfico de fração do espaço de planejamento para o projeto Box-Behnken.

A partir deste gráfico podemos ver que a variância relativa de um valor previsto é inferior a 0.35 em 50% da região de projeto. Este gráfico pode ser feito com a função FDSPlot que também está no pacote Vdgraph. Os gráficos de dispersão de variância são melhores para avaliar o quão próximo um projeto está de ser rotacionável e ter precisão uniforme.

Os gráficos de fração de espaço de planejamento são geralmente melhores para comparar projetos concorrentes para um problema de pesquisa.

library(Vdgraph)
FDSPlot(bbd3[ , 3:5], mod=2)
grid()

Figura 10.6: Gráfico de fração do espaço de planejamento para projeto Box-Behnken em três fatores.

Como exemplo de projeto Box-Behnken, considere o experimento realizado por Anderson (2003). Ele experimentou um modelo em escala de trabuco, um lançador de mísseis medieval, que foi distribuído a todos os estudantes de engenharia da South Dakota School of Mines and Technology para experimentação prática.

Ele queria determinar as configurações que lhe permitiriam lançar um objeto a uma determinada distância para se defender de invasores ou vendedores de porta em porta. Ele variou três pontos no trabuco. A, o comprimento do braço é de 4 a 8 polegadas da extremidade do contrapeso até o ponto onde os pesos foram pendurados; B, contrapeso de 10 a 20 libras; e C, peso do míssil de 2 a 3 onças.

O projeto Box-Behnken em níveis de fator reais é mostrado na Tabela 10.3 junto com a resposta, a distância que o míssil voou.

\[ \begin{array}{cccccccc}\hline \mbox{Corrida} & A & B & C & y \\\hline 1 & 4 & 10 & 2.5 & 33 \\ 2 & 8 & 10 & 2.5 & 85 \\ 3 & 4 & 20 & 2.5 & 86 \\ 4 & 8 & 20 & 2.5 & 113 \\ 5 & 4 & 15 & 2 & 75 \\ 6 & 8 & 15 & 2 & 105 \\ 7 & 4 & 15 & 3 & 40 \\ 8 & 8 & 15 & 3 & 89 \\ 9 & 6 & 10 & 2 & 83 \\ 10 & 6 & 20 & 2 & 108 \\ 11 & 6 & 10 & 3 & 49 \\ 12 & 6 & 20 & 3 & 101 \\ 13 & 6 & 15 & 2.5 & 88 \\ 14 & 6 & 15 & 2.5 & 91 \\ 15 & 6 & 15 & 2.5 & 91 \\\hline \end{array} \] Tabela 10.3: Resultados do projeto Box-Behnken para experimentos com trabuco.


10.3.3 Projeto composto pequeno


Quando o custo de cada experimento é alto e um experimentador está disposto a comprometer algumas das propriedades desejáveis de um projeto de superfície de resposta para reduzir o tamanho da execução, um pequeno projeto composto pode ser utilizado. Em um experimento composto pequeno, a parte \(2^k\) ou resolução V \(2^{k−p}\) do experimento composto central é substituída por uma resolução III \(2^{k−p}\).

No delineamento composto central, o \(2^k\) ou resolução V \(2^{k-p}\) foi incluído para permitir a estimatção de todos os termos de efeito principal linear e interações lineares por lineares de dois fatores. Hartley (1959) mostrou que quando pontos centrais e pontos axiais são incluídos, as interações de dois fatores podem ser estimadas com um fatorial fracionário de resolução III. Westlake (1965) obteve pequenos projetos compostos adicionais substituindo frações irregulares para a porção fatorial do experimento, e Draper (1985) e Draper and Lin (1990) substituíram os experimentos de Plackett-Burman pela porção fatorial para chegar a experimentos compostos ainda mais pequenos.

A Figura 10.7 é uma comparação gráfica do experimento composto central e composto pequeno para \(k = 2\) fatores. O experimento composto central tem cinco pontos centrais para um total de 13 execuções, enquanto o experimento composto pequeno tem um ponto central com um total de 7 execuções.

Figura 10.7: Comparação gráfica do planejamento composto central e do projeto composto pequeno, com \(I = AB\) para \(k=2\).

A Figura 10.8 é o gráfico de dispersão da variância para o experimento composto pequeno com \(k = 2\). Isso mostra que o experimento é quase rotacional dentro de um raio codificado de 0.8, mas não tem precisão uniforme, uma vez que a variância de previsão em escala é maior no centro do experimento do que em um raio codificado de 0.8.

O gráfico de dispersão da variância para o experimento composto central de precisão uniforme com \(k = 2\) seria muito semelhante à Figura 10.3.

library(Vdgraph)
data(SCDH2)
Vdgraph(SCDH2)
## number of design points= 7 
## number of factors= 2
##           Radius   Maximum  Minimum   Average
##  [1,] 0.00000000  7.000000 7.000000  7.000000
##  [2,] 0.07071068  6.970191 6.948635  6.959399
##  [3,] 0.14142136  6.881552 6.797610  6.839346
##  [4,] 0.21213203  6.736446 6.556125  6.645096
##  [5,] 0.28284271  6.538810 6.239520  6.385404
##  [6,] 0.35355339  6.294156 5.869268  6.072530
##  [7,] 0.42426407  6.009572 5.472977  5.722236
##  [8,] 0.49497475  5.693721 5.084391  5.353784
##  [9,] 0.56568542  5.356840 4.743389  4.989942
## [10,] 0.63639610  5.010741 4.474707  4.656979
## [11,] 0.70710678  4.668812 4.221746  4.384666
## [12,] 0.77781746  4.496707 3.987928  4.206277
## [13,] 0.84852814  4.867536 3.786901  4.158590
## [14,] 0.91923882  5.577372 3.633740  4.281883
## [15,] 0.98994949  6.702907 3.544785  4.619938
## [16,] 1.06066017  8.326965 3.537669  5.220039
## [17,] 1.13137085 10.538507 3.631328  6.132974
## [18,] 1.20208153 13.432630 3.846017  7.413032
## [19,] 1.27279221 17.110564 4.209751  9.118004
## [20,] 1.34350288 21.679676 4.767665 11.309186
## [21,] 1.41421356 27.253469 5.550248 14.051374
grid()

Figura 10.8: Gráfico de dispersão da variância de experimento composto pequeno com 2 fatores.

A Figura 10.9 é uma comparação dos gráficos de fração do espaço de projeto para os experimentos compostos centrais e compostos pequenos para \(k = 2\) fatores. Este gráfico é melhor para comparar os dois projetos. Pode-se observar, por referência à linha pontilhada horizontal superior no gráfico, que a variância relativa da predição seria inferior a 0.64 em 50% do espaço de projeto se o experimento composto pequeno fosse utilizado.

Por outro lado, a variância relativa da previsão seria inferior a 0.64 em mais de 80% da região do projeto se o experimento composto central fosse usado. A variância relativa da previsão é uniformemente menor para o experimento composto central. A escolha entre os dois projetos seria feita ponderando a precisão extra que o projeto composto central oferece em relação ao custo dos seis experimentos adicionais necessários.

# first get the ccd2
library(rsm)
ccd2 = ccd(2, n0 = c(3,2), alpha = "rotatable", randomize = FALSE)[ ,3:4]
Compare2FDS(SCDH2,ccd2,"SCD","CCD",mod=2)
grid()

Figura 10.9: Comparação de frações de parcelas de espaço de cálculo para o dimensionamento composto pequeno e o dimensionamento composto central para dois fatores


10.3.4 Projeto híbrido


Roquemore (1976) desenvolveu projetos híbridos que exigem ainda menos execuções do que os pequenos projetos compostos. Esses projetos foram construídos fazendo um projeto composto central em \(k-1\) fatores e adicionando um \(k\)-ésimo fator para que \(X^\top X\) tenha certas propriedades e o projeto seja quase rotacional. A Tabela 10.4 mostra um exemplo do projeto do Roquemore 310.

O experimento é rotulado como 310 porque tem 3 fatores e 10 execuções. Pode-se observar que existe um dimensionamento composto central nas colunas 1 e 2. Este dimensionamento é quase rotativo apenas dentro de um raio codificado de 0,5, conforme mostrado na Figura 10.10, mas como possui apenas 10 execuções, há zero graus de liberdade. para estimar o erro experimental.

\[ \begin{array}{cccccccc}\hline \mbox{Corrida} & x_1 & x_2 & x_3 \\\hline 1 & 0 & 0 & 1.2906 \\ 2 & 0 & 0 & -0.1360 \\ 3 & -1 & -1 & 0.6386 \\ 4 & 1 & -1 & 0.6386 \\ 5 & -1 & 1 & 0.6386 \\ 6 & 1 & 1 & 0.6386 \\ 7 & 1.7336 & 0 & -0.9273 \\ 8 & -1.7336 & 0 & -0.9273 \\ 9 & 0 & 1.7336 & -0.9273 \\ 10 & 0 & -1.736 & -0.9273 \\\hline \end{array} \] Tabela 10.4: Projeto Roquemore 310.

Alguns dos projetos híbridos de Roquemore deixam um grau de liberdade para estimar o erro experimental, mas nenhum deles possui réplicas puras que são necessárias para testar a adequação do modelo quadrático geral. Portanto, esses projetos só devem ser usados quando o experimentador estiver confiante de que o modelo quadrático geral representará bem os dados e houver uma estimativa independente do erro experimental a partir de dados anteriores.

library(Vdgraph)
data(D310)
Vdgraph(D310)
## number of design points= 10 
## number of factors= 3
##           Radius   Maximum   Minimum   Average
##  [1,] 0.00000000  9.334145  9.334145  9.334145
##  [2,] 0.08660254  9.836354  8.748894  9.288908
##  [3,] 0.17320508 10.263005  8.086022  9.155941
##  [4,] 0.25980762 10.634480  7.363846  8.943485
##  [5,] 0.34641016 10.984065  6.613588  8.665267
##  [6,] 0.43301270 11.357945  5.879370  8.340511
##  [7,] 0.51961524 11.815205  5.218214  7.993929
##  [8,] 0.60621778 12.427835  4.700045  7.655728
##  [9,] 0.69282032 13.280723  4.407687  7.361606
## [10,] 0.77942286 14.471657  4.316738  7.152753
## [11,] 0.86602540 16.111330  4.336948  7.075849
## [12,] 0.95262794 18.323332  4.508247  7.183069
## [13,] 1.03923048 21.244157  4.875242  7.532078
## [14,] 1.12583302 25.023198  5.486722  8.186035
## [15,] 1.21243557 29.822749  5.751032  9.213588
## [16,] 1.29903811 35.818008  5.798200 10.688880
## [17,] 1.38564065 43.197071  6.047784 12.691543
## [18,] 1.47224319 52.160936  6.554966 15.306704
## [19,] 1.55884573 62.923501  7.377301 18.624980
## [20,] 1.64544827 75.711567  8.575243 22.742480
## [21,] 1.73205081 90.764835 10.212732 27.760806
grid()

Figura 10.10: Gráfico de dispersão da variância para o projeto Roquemore 310.


10.4 Criando planejamentos de superfície de resposta padrão em R


Projetos de superfície de resposta padrão podem ser criados usando pacotes R. O pacote rsm possui funções para criar projetos compostos centrais e projetos Box-Behnken. Este pacote também tem uma função para listar vários planejamentos compostos centrais disponíveis para um número específico de fatores.

O pacote Vdgraph contém conjuntos de dados armazenados que contêm pequenos projetos compostos e híbridos para 2 a 6 fatores.

Se um projeto de superfície de resposta padrão para três fatores for desejado, o código R abaixo ilustra a chamada à função ccd.pick no pacote rsm para mostrar os possíveis projetos compostos centrais disponíveis.

library(rsm)
ccd.pick(k=3)
##    n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s  N alpha.rot alpha.orth
## 1    8    9      1   6    6     1     1     1 29  1.681793   1.680336
## 2    8    2      1   6    1     1     1     1 17  1.681793   1.673320
## 3    8    6      1   6    4     1     1     1 24  1.681793   1.690309
## 4    8    5      1   6    3     1     1     1 22  1.681793   1.664101
## 5    8   10      1   6    7     1     1     1 31  1.681793   1.699673
## 6    8    8      1   6    5     1     1     1 27  1.681793   1.658312
## 7    8    3      1   6    2     1     1     1 19  1.681793   1.705606
## 8    8    7      1   6    5     1     1     1 26  1.681793   1.712698
## 9    8    4      1   6    2     1     1     1 20  1.681793   1.632993
## 10   8    4      1   6    3     1     1     1 21  1.681793   1.732051

Na saída acima, 10 projetos possíveis são mostrados. Nesta saída, n.c representa o número de pontos na porção fatorial do experimento, neste caso um fatorial completo \(2^3\), n0.c representa o número de pontos centrais na porção fatorial do experimento, blks.c representa o número de blocos no fatorial parte do projeto, n.s representa o número de pontos na porção axial do projeto, n0.s representa o número de pontos centrais na porção axial do projeto, bbr.c, wbr.s e bbr.s representam a número de cópias da porção fatorial, o número de réplicas de cada ponto axial e o número de cópias da porção axial do plano, respectivamente.

Finalmente, N representa o número total de execuções no projeto, alpha.rot e alpha.orth representam o valor em unidades codificadas para um projeto rotacional ou ortogonal.

Normalmente, a rotabilidade deve ser escolhida se um compósito central for usado e o número de pontos centrais pode ser ajustado para fornecer uma precisão quase uniforme. Quando projetos compostos centrais são bloqueados, a ortogonalidade deve ser escolhida para tornar os efeitos de bloco ortogonais ao modelo.

Para criar um dos planos listados na saída da função ccd.pick, a função ccd no pacote rsm pode ser usada conforme mostrado abaixo.

library(rsm)
ccd.up = ccd(y~x1+x2+x3,n0=c(4,2),alpha="rotatable", randomize=FALSE)
head(ccd.up)
##   run.order std.order x1.as.is x2.as.is x3.as.is  y Block
## 1         1         1       -1       -1       -1 NA     1
## 2         2         2        1       -1       -1 NA     1
## 3         3         3       -1        1       -1 NA     1
## 4         4         4        1        1       -1 NA     1
## 5         5         5       -1       -1        1 NA     1
## 6         6         6        1       -1        1 NA     1
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is


Este código produz o nono plano que possui n0.c=4 e n0.s=2. Este projeto possui precisão uniforme como pode ser visualizado pelo seu gráfico de dispersão da variância produzido com o comando abaixo.

library(Vdgraph)
Vdgraph(ccd.up[ , 3:5])
## number of design points= 20 
## number of factors= 3
##           Radius   Maximum   Minimum   Average
##  [1,] 0.00000000  3.326805  3.326805  3.326805
##  [2,] 0.08660254  3.320828  3.320828  3.320828
##  [3,] 0.17320508  3.303837  3.303837  3.303837
##  [4,] 0.25980762  3.278640  3.278640  3.278640
##  [5,] 0.34641016  3.249923  3.249923  3.249923
##  [6,] 0.43301270  3.224241  3.224241  3.224241
##  [7,] 0.51961524  3.210026  3.210026  3.210026
##  [8,] 0.60621778  3.217583  3.217583  3.217583
##  [9,] 0.69282032  3.259089  3.259089  3.259089
## [10,] 0.77942286  3.348596  3.348596  3.348596
## [11,] 0.86602540  3.502029  3.502029  3.502029
## [12,] 0.95262794  3.737186  3.737186  3.737186
## [13,] 1.03923048  4.073740  4.073740  4.073740
## [14,] 1.12583302  4.533236  4.533236  4.533236
## [15,] 1.21243557  5.139093  5.139093  5.139093
## [16,] 1.29903811  5.916603  5.916603  5.916603
## [17,] 1.38564065  6.892934  6.892934  6.892934
## [18,] 1.47224319  8.097125  8.097125  8.097125
## [19,] 1.55884573  9.560089  9.560089  9.560089
## [20,] 1.64544827 11.314613 11.314613 11.314613
## [21,] 1.73205081 13.395357 13.395357 13.395357
grid()


O plano ccd.up produzido no código acima possui níveis de fator codificados que variam entre -1.68 e +1.68. Se a opção de codificação for utilizada, conforme mostrado no código abaixo, a listagem resultante do experimento mostra os níveis reais dos fatores.

ccd.up = ccd(y~x1+x2+x3,n0=c(4,2),alpha="rotatable", 
             coding=list(x1~(Temp-150)/10,x2~(Press-50)/5, x3~(Rate-4/1)), randomize=FALSE)
head(ccd.up)
##   run.order std.order Temp Press Rate  y Block
## 1         1         1  140    45    3 NA     1
## 2         2         2  160    45    3 NA     1
## 3         3         3  140    55    3 NA     1
## 4         4         4  160    55    3 NA     1
## 5         5         5  140    45    5 NA     1
## 6         6         6  160    45    5 NA     1
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Temp - 150)/10
## x2 ~ (Press - 50)/5
## x3 ~ (Rate - 4/1)


No entanto, eles ainda são armazenados em níveis codificados para fins computacionais ao ajustar um modelo ao projeto e examinar a superfície ajustada. Outras opções para alpha são “spherical” para colocar pontos axiais na mesma distância da origem que os cantos do desenho e “faces” para colocar pontos axiais em níveis fatoriais altos e baixos, às vezes chamado de desenho de cubo centrado na face.

A função bbd no pacote rsm gera projetos Box-Behnken em três a sete fatores. Por exemplo, o código abaixo produz o projeto mostrado na Tabela 10.3.

library(rsm)
Treb = bbd(y~x1+x2+x3, randomize=FALSE, n0=3, 
           coding=list(x1~(A-6)/2, x2~(B-15)/5, x3~(C-2.5)/.5))
Treb
##    run.order std.order A  B   C  y
## 1          1         1 4 10 2.5 NA
## 2          2         2 8 10 2.5 NA
## 3          3         3 4 20 2.5 NA
## 4          4         4 8 20 2.5 NA
## 5          5         5 4 15 2.0 NA
## 6          6         6 8 15 2.0 NA
## 7          7         7 4 15 3.0 NA
## 8          8         8 8 15 3.0 NA
## 9          9         9 6 10 2.0 NA
## 10        10        10 6 20 2.0 NA
## 11        11        11 6 10 3.0 NA
## 12        12        12 6 20 3.0 NA
## 13        13        13 6 15 2.5 NA
## 14        14        14 6 15 2.5 NA
## 15        15        15 6 15 2.5 NA
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (A - 6)/2
## x2 ~ (B - 15)/5
## x3 ~ (C - 2.5)/0.5


O pacote rsm não possui funções para criar projetos híbridos ou pequenos projetos compostos diretamente, mas o pacote Vdgraph possui os projetos mais comuns armazenados como conjuntos ou quadros de dados que podem ser recuperados.

A Tabela 10.5, mostrada a continuação, lista os planos disponíveis como conjuntos de dados neste pacote.

\[ \begin{array}{ll} \mbox{Projeto composto pequeno} & \\ \mbox{Nome do conjunto de dados} & \mbox{Descrição} \\ \hline \mbox{SCDDL5} & \mbox{Projeto para 5 fatores de Draper and Lin} \\ \mbox{SCDH2} & \mbox{Projeto para dois fatores de Hartley} \\ \mbox{SCDH3} & \mbox{Projeto para treis fatores de Hartley} \\ \mbox{SCDH4} & \mbox{Projeto para quatro fatores de Hartley} \\ \mbox{SCDH5} & \mbox{Projeto para cinco fatores de Hartley} \\ \mbox{SCDH62} & \mbox{Projeto para seis fatores de Hartley} \\ & & \\ \mbox{Projeto híbrido} & \\ \mbox{Nome do conjunto de dados} & \mbox{Descrição} \\ \hline \mbox{D310} & \mbox{Projeto híbrido de Roquemore D310} \\ \mbox{D311A} & \mbox{Projeto híbrido de Roquemore D311A} \\ \mbox{D311B} & \mbox{Projeto híbrido de Roquemore DD311B} \\ \mbox{D416A} & \mbox{Projeto híbrido de Roquemore D416A} \\ \mbox{D416B} & \mbox{Projeto híbrido de Roquemore D416B} \\ \mbox{D416C} & \mbox{Projeto híbrido de Roquemore D416C} \\ \mbox{D628A} & \mbox{Projeto híbrido de Roquemore D628A} \\ & & \\ \mbox{Projeto hexagonal} & \\ \mbox{Nome do conjunto de dados} & \mbox{Descrição} \\ \hline \mbox{Hex2} & \mbox{Projeto hexagonal em dois fatores} \\ \end{array} \] Tabela 10.5: Projetos de superfície de resposta de execução reduzida no pacote Vdgraph.

O Hex2 é um projeto rotacional de precisão quase uniforme em dois fatores que contém apenas nove execuções. Essas corridas estão situadas nos vértices de um hexágono em um espaço bidimensional.

Um projeto pode ser recuperado conforme mostrado no exemplo abaixo.

library(Vdgraph)
data(D310)
D310
##         x1     x2      x3
## 1   0.0000  0.000  1.2906
## 2   0.0000  0.000 -0.1360
## 3  -1.0000 -1.000  0.6386
## 4   1.0000 -1.000  0.6386
## 5  -1.0000  1.000  0.6386
## 6   1.0000  1.000  0.6386
## 7   1.7636  0.000 -0.9273
## 8  -1.7636  0.000 -0.9273
## 9   0.0000  1.736 -0.9273
## 10  0.0000 -1.736 -0.9273


Os conjuntos de dados no Vdgraph possuem os fatores denominados \(x_1\) - \(x_k\) e níveis de fator codificados. Para fazer uma listagem do experimento em ordem aleatória e nos níveis reais dos fatores, use os comandos semelhantes aos mostrados abaixo.

des = transform(D310, Temp=10*x1+150, Press=5*x2+50, Rate=x3+4)
des[sample(1:10) ,4:6]
##       Temp Press   Rate
## 9  150.000 58.68 3.0727
## 3  140.000 45.00 4.6386
## 6  160.000 55.00 4.6386
## 1  150.000 50.00 5.2906
## 5  140.000 55.00 4.6386
## 2  150.000 50.00 3.8640
## 7  167.636 50.00 3.0727
## 10 150.000 41.32 3.0727
## 8  132.364 50.00 3.0727
## 4  160.000 45.00 4.6386

10.5 Projetos não padrões de superfície de resposta


Algumas situações de projeto não se prestam ao uso de projetos de superfície de resposta padrão por razões como

  1. a região de experimentação tem formato irregular,

  2. nem todas as combinações dos níveis de fator são viáveis,

  3. não há um modeo não-linear ou um modelo não-padrão linear.

Nessas situações, os projetos de superfície de resposta padrão não serão apropriados.

Uma maneira de construir um projeto de superfície de resposta nesses casos é usar um algoritmo de computador para construir um projeto I-ótimo.

Para um exemplo da primeira situação em que a região de experimentação tem formato irregular, Atkinson et al. (2007) descrevem um experimento para investigar o desempenho de um motor de combustão interna. Dois fatores em estudo foram o avanço da faísca e o torque. Ambos são independentemente variáveis, portanto a região de desenho codificada é um quadrado.

Contudo, para combinações de fatores fora da região pentagonal mostrada na Figura 10.11, é provável que o motor não funcione ou seja seriamente danificado. Portanto, a região experimental é irregular e os projetos padrão de superfície de resposta de dois fatores não funcionariam nesta região.

!{}(Figura1011.png) Figura 10.11: Região experimental para experimento de motor.

Um projeto de superfície de resposta pode ser construído para este problema definindo pontos candidatos que se enquadram em uma grade sobre a região experimental, mostrada na Figura 10.11. Isso pode ser feito no R e a função optFederov no pacote AlgDesign pode ser usada para selecionar um subconjunto I-ótimo para ajustar o modelo quadrático geral.

Como exemplo da segunda situação em que nem todas as combinações possíveis do factor são viáveis, considere um estudo típico em QSAR (Quantitative Structure Activity Relation) utilizado na concepção de medicamentos. Neste tipo de estudo, quando é descoberto um composto líder com bioatividade desejável, muitos derivados deste composto são considerados na tentativa de encontrar um que aumente o efeito desejado.

Por exemplo, a Figura 10.12 mostra a estrutura geral das hidroxifenilureias, que demonstraram ser ativas como antioxidantes.

Figura 10.12: Estrutura geral das hidroxifenilureias.

O composto possui quatro locais receptores R, R’, R’’ e R’’’, mostrados na Figura 10.12, onde diferentes átomos ou moléculas substituintes podem ser ligadas para modificar a estrutura básica.

Com muitas moléculas substituintes possíveis que podem ser adicionadas em cada sítio receptor, existe uma grande biblioteca de candidatos que podem ser sintetizados e testados quanto à atividade. A Tabela 10.6 é uma lista ilustrativa que mostra 36 possíveis hidroxifenilureias diferentes de uma lista mais longa considerada por Deeb et al. (2008).

Esta lista de candidatos é armazenada em um conjunto de dados no pacote R daewr e pode ser recuperada com o comando abaixo.

library(daewr)
data(qsar)


Para cada composto variante possível na biblioteca existem vários descritores físico-químicos que podem ser calculados. Na Tabela 10.6 são mostrados três descritores: Energia de Hidratação (HE), Momento de Dipolo Molecular na direção z (DMz) e Índice de Simetria (S0K).

\[ \begin{array}{cccccccc}\hline \mbox{Composto} & \mbox{R} & \mbox{R'} & \mbox{R''} & \mbox{R'''} & \mbox{HE} & \mbox{DMz} & \mbox{S0K} \\\hline 1 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{CH}_3 & -12.221 & -0.162 & 64.138 \\ 2 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -14.015 & -0.068 & 88.547 \\ 3 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{Ph} & -14.502 & -0.372 & 85.567 \\ 4 & \mbox{H} & \mbox{H} & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -14.893 & 1.035 & 96.053 \\ 5 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_3 & -12.855 & 1.091 & 74.124 \\ 6 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -14.628 & 1.115 & 99.002 \\ 7 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{Ph} & -15.123 & 1.554 & 96.053 \\ 8 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -15.492 & 2.221 & 106.607 \\ 9 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{CH}_3 & -11.813 & 1.219 & 77.020 \\ 10 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -13.593 & 1.188 & 101.978 \\ 11 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{Ph} & -14.088 & 1.621 & 99.002 \\ 12 & \mbox{CH}_3 & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -14.460 & 2.266 & 109.535 \\ 13 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{CH}_3 & -8.519 & -0.560 & 71.949 \\ 14 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{CH}_2\mbox{Ph} & -10.287 & -0.675 & 96.600 \\ 15 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{Ph} & -10.798 & -0.134 & 96.620 \\ 16 & \mbox{CH}_3 & \mbox{H}_3 & \mbox{CH}_3 & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -11.167 & 0.418 & 104.047 \\ 17 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{CH}_3 & -12.245 & -0.609 & 67.054 \\ 18 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -13.980 & -0.561 & 88.547 \\ 19 & \mbox{H} & \mbox{H} & \mbox{H} & \mbox{Ph} & -14.491 & -0.561 & 88.547 \\ 20 & \mbox{H} & \mbox{H} & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -14.888 & 1.478 & 99.002 \\ 21 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_3 & -11.414 & -1.888 & 77.020 \\ 22 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -13.121 & -1.692 & 101.978 \\ 23 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & \mbox{Ph} & -13.660 & -1.893 & 99.002 \\ 24 & \mbox{H} & \mbox{OCH}_3 & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -14.012 & -2.714 & 109.535 \\ 25 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{CH}_3 & -10.029 & -1.891 & 79.942 \\ 26 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -11.740 & -1.652 & 104.977 \\ 27 & \mbox{H} & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & \mbox{Ph} & -12.329 & -1.902 & 101.978 \\ 28 & \mbox{OCH}_3 & \mbox{OC}_2\mbox{H}_5 & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -12.637 & -2.762 & 112.492 \\ 29 & \mbox{OCH}_3 & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_3 & -12.118 & -2.994 & 81.106 \\ 30 & \mbox{OCH}_3 & \mbox{OCH}_3 & \mbox{H} & \mbox{CH}_2\mbox{Ph} & -13.892 & -2.845 & 106.299 \\ 31 & \mbox{OCH}_3 & \mbox{OCH}_3 & \mbox{H} & \mbox{Ph} & -14.456 & -2.926 & 103.230 \\ 32 & \mbox{OCH}_3 & \mbox{OCH}_3 & \mbox{H} & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -14.804 & -3.780 & 113.856 \\ 33 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{CH}_3 & -9.209 & -0.423 & 74.871 \\ 34 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{CH}_2\mbox{Ph} & -10.970 & -0.302 & 99.603 \\ 35 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & \mbox{Ph} & -11.488 & -0.453 & 96.600 \\ 36 & \mbox{CH}_3 & \mbox{H} & \mbox{CH}_3 & 2\mbox{CH}_3\mbox{OC}_6\mbox{H}_4 & -11.868 & -1.322 & 107.010 \\\hline \end{array} \] Tabela 10.6: Biblioteca de compostos de hidroxifenilureia substituintes.

O objetivo de um estudo QSAR seria

  1. sintetizar um subconjunto de moléculas da grande biblioteca e testá-las para determinar sua capacidade de eliminar FRs de oxigênio, conforme medido pela constante de ligação \(\log\big(K_{app}\big)\), e

  2. um modelo relacionando a resposta \(\log\big(K_{app}\big)\) para vários descritores físico-químicos, neste caso, as variáveis HE, DMz e S0K.

Uma vez encontrado um modelo que relacione a resposta em função das variáveis independentes, a combinação das variáveis independentes que se prevê fornecer a resposta máxima pode ser determinada. Os compostos na grande biblioteca que possuem valores dos descritores físico-químicos (variáveis independentes) mais próximos daqueles previstos para resultar em uma resposta máxima podem então ser sintetizados e testados quanto à atividade. Desta forma, são descobertos compostos variantes que têm actividade aumentada.

Nos estudos QSAR, os planejamentos de superfícies de resposta padrão não podem ser utilizados porque nem todas as combinações potenciais das variáveis independentes são possíveis. Apenas um subconjunto das combinações de variáveis independentes que existem na grande biblioteca pode ser escolhido em um subconjunto ou plano experimental a ser testado.

Novamente, a função optFederov no pacote AlgDesign pode ser utilizada para criar um subconjunto I-ótimo ou D-ótimo da biblioteca maior. Neste caso, a lista de candidatos será composta por todos os compostos da biblioteca. O plano deve consistir em um subconjunto dos compostos da biblioteca.

Se o modelo quadrático geral for usado para testar os dados com três fatores, existem 10 coeficientes no modelo quadrático geral, portanto o número de execuções no projeto deve ser pelo menos \(n = 10\). No código R abaixo da opção nRepeats= 40 instrui a função optFederov a fazer 40 pesquisas diferentes para o experimento D-ótimo com um experimento inicial aleatório e manter os melhores resultados.

A opção nTrials=16 instrui a função optFederov a fazer um experimento com 16 execuções. O plano será armazenado no arquivo de dados (data frame) desgn1$design, o quarto item da lista desgn1. A chamada da função optFederov é repetida alterando a opção de critério de D para I para criar um plano I-ótimo que é armazenado no conjunto de dados desgn2$design.

A última instrução chama a função Compare2FDS para comparar a fração dos gráficos do espaço de design para os dois experimentos.

library(daewr)
data(qsar)
library(AlgDesign)
desgn1 = optFederov( ~quad(.), data=qsar, nTrials=15, center=TRUE, criterion="D", nRepeats=40)
desgn2 = optFederov( ~quad(.), data=qsar, nTrials=15, center=TRUE, criterion="I", nRepeats=40)
library(Vdgraph)
Compare2FDS(desgn1$design, desgn2$design, "D-optimal", "I-optimal", mod=2)
grid()

Figura 10.13: Comparação dos gráficos de fração de espaço do projeto D-Ótimo e o projeto I-ótimo.

desgn2$design
##    Compound      HE    DMz     S0K
## 1         1 -12.221 -0.162  64.138
## 4         4 -14.893  1.035  96.053
## 9         9 -11.813  1.219  77.020
## 12       12 -14.460  2.266 109.535
## 13       13  -8.519 -0.560  71.949
## 14       14 -10.287 -0.675  96.600
## 16       16 -11.167  0.418 104.047
## 19       19 -14.491 -0.561  88.547
## 22       22 -13.121 -1.692 101.978
## 28       28 -12.637 -2.762 112.492
## 29       29 -12.118 -2.994  81.106
## 32       32 -14.804 -3.780 113.856
## 33       33  -9.209 -0.423  74.871
## 34       34 -10.970 -0.302  99.603
## 36       36 -11.868 -1.322 107.010


O plano I-ótimo pode ser listado conforme mostrado acima. Depois que os projetos são criados com a função optFederov, eles podem ser comparados em relação à variação de um valor previsto. Embora nenhum dos projetos seja rotacional, o projeto I-ótimo é ligeiramente melhor, uma vez que sua variância relativa de predição aparece uniformemente menor do que a variância relativa de predição do projeto D-ótimo, como mostrado na Figura 10.13.

A terceira situação em que um projeto de superfície de resposta não padronizado deve ser usado é quando o modelo não é padronizado ou não linear. Considere modelar o metabolismo da tetraciclina, um antibiótico comum. A Figura 10.14 mostra um diagrama do modelo de dois compartimentos útil em farmacocinética para representar o metabolismo de medicamentos.

Figura 10.14: Diagrama do modelo de dois compartimentos para metabolismo da tetraciclina.

Quando o medicamento é administrado por via oral:

  1. \(\gamma_0\) representa a concentração inicial no trato gastrointestinal (GI),

  2. \(\gamma_1(x)\) representa a concentração no sangue no tempo \(x\),

  3. \(\kappa_1\) representa a velocidade constante na qual o medicamento se move do trato GI para o sangue,

  4. \(\kappa_2\) é a velocidade constante na qual o medicamento é eliminado do sangue.

Com base na cinética química, uma equação relacionando a concentração no sangue no tempo \(x\) pode ser derivada e é mostrada na Equação (10.7) \[\begin{equation*} \tag{10.7} y=\gamma_1(x)=\gamma_0\Big(e^{-\kappa_1(x-t_0)}-e^{-\kappa_2(x-t_0)} \Big), \end{equation*}\] onde \(t_0\) é o tempo inicia (dead time).

Este é um exemplo de modelo mecanicista derivado de princípios físicos. Para determinar os parâmetros no modelo \((\gamma_0, \kappa_1, \kappa_2, t_0)\), um sujeito toma uma dose oral do medicamento e pequenas amostras de sangue são coletadas nos tempos \(x_i\), \(i = 1,\cdots,n\).

A resposta ou concentração do medicamento no sangue \(\gamma_1(x_i)\) é determinada em cada tempo de amostragem ou nível de fator e o modelo 10.7 é ajustado aos dados. Uma vez ajustado, o modelo mecanístico pode ser usado para prever valores da resposta fora da faixa experimental de \(x_i\), \(i = 1,\cdots,n\) já que é baseado em princípios físicos. Isto não é verdade para o modelo quadrático geral (10.4), uma vez que é apenas uma aproximação por Série de Taylor da equação verdadeira em torno de algum ponto \((x)\).

Neste problema existe apenas uma variável independente e uma resposta e existem quatro parâmetros no modelo. Como o modelo é não linear nos parâmetros, um projeto simples de quatro tempos igualmente espaçados para coletar amostras de sangue, não seria eficiente para estimar os parâmetros. Em vez disso, uma pesquisa I-ótima pode ser feita para determinar os pontos de dimensionamento.

O primeiro passo na criação de um projeto I-ótimo para um modelo não linear é linearizar o modelo em torno de algum ponto, conforme descrito por Bates and Watts (2007).

Por exemplo, linearizando \(f(x;\gamma_0,\kappa_1,\kappa_2,t_0)\) acerca do ponto \((\gamma_0^*,\kappa_1^*,\kappa_2^*,t_0^*)\) \[\begin{equation*} \tag{10.8} f(x;\gamma_0,\kappa_1,\kappa_2,t_0) \approx f(x;\gamma_0^*,\kappa_1^*,\kappa_2^*,t_0^*)+(\gamma_0-\gamma_0^*)\left.\dfrac{\partial f}{\partial\gamma_0}\right|_{\gamma_0=\gamma_0^*}+(\kappa_1-\kappa_1^*)\left.\dfrac{\partial f}{\partial\kappa_1}\right|_{\kappa_a=\kappa_1^*} \\ \qquad \qquad \qquad \qquad \quad +(\kappa_2-\kappa_2^*)\left.\dfrac{\partial f}{\partial\kappa_2}\right|_{\kappa_2=\kappa_2^*} +(t_0-t_0^*)\left.\dfrac{\partial f}{\partial t_0}\right|_{t_0=t_0^*}, \end{equation*}\] que é linear nas variáveis \(\partial f/\partial \gamma_0\), \(\partial f/\partial \kappa_1\), \(\partial f/\partial \kappa_2\) e \(\partial f/\partial t_0\), que são todas funções de \(x\).

Para o modelo de compartimento na Equação (10.7) \[ \dfrac{\partial f}{\partial\gamma_0} = e^{-\kappa_1(x-t_0)}-e^{-\kappa_2(x-t_0)}, \\ \dfrac{\partial f}{\partial\kappa_1} = -\gamma_0(x-t_0)e^{-\kappa_1(x-t_0)}, \\ \dfrac{\partial f}{\partial\gamma_0} = -\gamma_0(x-t_0)e^{-\kappa_2(x-t_0)}, \\ \dfrac{\partial f}{\partial\gamma_0} = \gamma_0\kappa_1 e^{-\kappa_1(x-t_0)}-\gamma_0\kappa_2 e^{-\kappa_2(x-t_0)}\cdot \]

A estratégia é criar uma grade de candidatos na variável independente \(x\), calcular os valores das quatro derivadas parciais usando estimativas iniciais dos valores dos parâmetros em cada ponto candidato e, em seguida, usar a função optFederov para selecionar um subconjunto I-ótimo da grade.

O código R abaixo cria a grade e avalia as derivadas parciais usando as estimativas iniciais dos parâmetros \[ \gamma_0^0 =2.65, \kappa_1^0 =0.15, \kappa_2^0 = 0.72 \quad \mbox{e} \quad t_0^0 =0.41\cdot \]

k1 = 0.15; k2 = 0.72; gamma0 = 2.65; t0 = 0.41
x = c(seq(1:25))
dfdk1 = c(rep(0, 25))
dfdk2 = c(rep(0, 25))
dfdgamma0 = c(rep(0, 25))
dfdt0 = c(rep(0, 25))
for (i in 1:25) { 
  dfdk1[i] = -1 * gamma0 * exp(-k1 * (x[i] - t0)) *(x[i] - t0) 
  dfdk2[i] = gamma0 * exp(-k2 * (x[i] - t0)) * (x[i] - t0)
  dfdgamma0[i] = exp(-k1 * (x[i] - t0)) - exp( -k2 * ( x[i] - t0))
  dfdt0[i] = gamma0 * exp(-k1 * (x[i] - t0)) * k1 - gamma0 * exp(-k2 * (x[i] - t0)) * k2; }
grid = data.frame(x, dfdk1, dfdk2, dfdgamma0, dfdt0)
library(AlgDesign)
desgn2 = optFederov( ~ -1 + dfdk1 + dfdk2 + dfdgamma0 + dfdt0, data = grid, 
                     nTrials = 4, center = TRUE, criterion = "I", nRepeats = 20)
desgn2$design
##     x     dfdk1        dfdk2  dfdgamma0        dfdt0
## 1   1 -1.431076 1.022374e+00 0.26140256 -0.883809267
## 2   2 -3.319432 1.341105e+00 0.46952112 -0.294138728
## 5   5 -6.110079 4.464802e-01 0.46562245  0.129639675
## 25 25 -1.629706 1.333237e-06 0.02500947  0.009941233


No código acima, a função optFederov é usada para criar um planejamento que consiste em um subconjunto da grade. Como existem \(p = 4\) parâmetros no modelo, apenas quatro pontos do plano serão necessários para estimar os parâmetros, portanto a opção nTrials = 4 é usada na chamada à função optFederov. Um experimentador pode então replicar cada ponto do projeto quantas vezes forem necessárias para obter uma estimativa do erro experimental.

O código acima resulta em quatro tempos distintos para a coleta de amostras de sangue, 1, 2, 5 e 25 e, devido à forma do modelo e às estimativas iniciais dos parâmetros, esses pontos estão longe de ser igualmente espaçados. Como a função optFederov nem sempre pode encontrar o plano ideal e do fato de que os planos ideais podem não ser únicos, executar esse código novamente pode resultar em um plano ligeiramente diferente.

Conforme afirmado por Bates and Watts (2007), a eficiência do projeto do modelo não linear dependerá do estágio em que o pesquisador se encontra na investigação. Se a forma da equação for conhecida, mas não os valores dos parâmetros, o projeto criado pela função optFederov será chamado de projeto inicial.

Depois que os dados forem coletados de um projeto inicial e o modelo for ajustado, as estimativas dos parâmetros serão uma melhoria nas estimativas iniciais. Portanto, as estimativas iniciais dos valores dos parâmetros podem ser substituídas pelas estimativas e o projeto inicial pode ser aumentado com \(p = 4\) pontos de projeto adicionais usando a opção augment = TRUE na chamada de função semelhante ao exemplo mostrado na Seção 6.5.2.


10.6 Ajustando o modelo de superfície de resposta no R


Esta seção descreverá como ajustar um modelo de superfície de resposta usando o pacote R rsm e um modelo mecanístico não linear usando a função R nls.


10.6.1 Ajustando um modelo linear e verificando a curvatura


Se a porção fatorial de um planejamento composto central juntamente com os pontos centrais for concluída em um bloco inicial de experimentos conforme descrito na Seção 10.3.1, o modelo linear (10.6) deverá ser ajustado aos dados e verificado quanto à adequação. Se o modelo linear for adequado, os pontos axiais e os pontos centrais adicionais no dimensionamento composto central não serão necessários.

Para verificar a adequação, as somas residuais dos quadrados do modelo 10.6 devem ser particionadas na porção devido à replicação pura e na porção devido ao desvio quadrático do modelo. Isso pode ser feito usando a função rsm no pacote R homônimo. Por exemplo, o conjunto de dados cement no pacote daewr contém os dados da Tabela 10.1.

Este conjunto de dados contém uma variável adicional Block não mostrada na Tabela 10.1. Block = 1 para as primeiras 11 execuções, fatorial + pontos centrais, na Tabela 10.1 e Block = 2 para as últimas nove execuções, axiais + pontos centrais.

No código R abaixo, os dados são recuperados do pacote daewr. Esse conjuto de dados foi criado com a função ccd no pacote rsm, portanto, os fatores são armazenados internamente como os fatores codificados \(x_1\), \(x_2\) e \(x_3\).

library(daewr)
data(cement)
head(cement)
##      Block WatCem BlackL  SNF     y
## C1.1     1   0.33   0.12 0.08 109.5
## C1.2     1   0.35   0.12 0.08 117.0
## C1.3     1   0.33   0.18 0.08 110.5
## C1.4     1   0.35   0.18 0.08 121.0
## C1.5     1   0.33   0.12 0.12 120.0
## C1.6     1   0.35   0.12 0.12 130.0
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (WatCem - 0.34)/0.01
## x2 ~ (BlackL - 0.15)/0.03
## x3 ~ (SNF - 0.1)/0.02


Uma análise é realizada nos dados do bloco 1 usando R. A fórmula SO(x1,x2,x3) na chamada para a função rsm cria o modelo quadrático completo (10.4). No entanto, quando os pontos axiais são deixados de fora do planejamento composto central, como no caso do bloco 1, todos os termos quadráticos ficam confusos.

O teste no termo PQ(x1,x2,x3) na tabela ANOVA resultante é um teste dos termos quadráticos confundidos.

library(rsm)
grout.lin = rsm(y ~ SO(x1, x2, x3),data = cement, subset = (Block == 1))
anova(grout.lin)
## Analysis of Variance Table
## 
## Response: y
##                 Df Sum Sq Mean Sq F value   Pr(>F)   
## FO(x1, x2, x3)   3 465.13 155.042 80.3094 0.002307 **
## TWI(x1, x2, x3)  3   0.25   0.083  0.0432 0.985889   
## PQ(x1, x2, x3)   1  37.88  37.879 19.6207 0.021377 * 
## Residuals        3   5.79   1.931                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


As somas dos quadrados para a parte linear, rotulado FO(x1, x2, x3) na saída, leva em conta simultaneamente os três termos lineares no modelo. As somas dos quadrados do produto vetorial, rotulado TWI(x1, x2, x3) na saída, leva em conta simultaneamente os três termos de interação no modelo.

As somas dos quadrados para a parte quadrática, rotulado PQ(x1, x2, x3) na saída, leva em conta os termos quadráticos confusos ou o desvio da linearidade. Finalmente, as somas dos quadrados Residual representam as somas puras dos erros dos quadrados devido aos pontos centrais replicados.

Ao ajustar o modelo linear com interações (10.6) com a função R lm, as somas dos quadrados para parte quadrática seriam combinadas com as somas dos quadrados dos erros, resultando em quatro graus de liberdade para o erro.

O teste \(F\) no termo quadrático é um teste de adequação do modelo linear. Como é significativo no nível \(\alpha = 0.05\) na tabela acima, indica que o desvio quadrático do modelo linear é significativo e, portanto, o modelo linear (10.6) não é adequado para representar os primeiros 11 pontos de dados em Tabela 10.1.

Portanto, se os experimentos mostrados na Tabela 10.1 fossem executados em dois blocos, seria necessário executar o segundo bloco que inclui os pontos axiais.


10.6.2 Ajustando o modelo quadrático geral


A função rsm é mais útil para ajustar o modelo quadrático geral e analisar a superfície ajustada. Por exemplo, considere ajustar o modelo quadrático geral aos dados do trabuco na Tabela 10.3.

Os comandos para recuperar o conjunto de dados do pacote daewr são:

library(daewr)
data(Treb)
head(Treb)
##   A  B   C   y
## 1 4 10 2.5  33
## 2 8 10 2.5  85
## 3 4 20 2.5  86
## 4 8 20 2.5 113
## 5 4 15 2.0  75
## 6 8 15 2.0 105
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (A - 6)/2
## x2 ~ (B - 15)/5
## x3 ~ (C - 2.5)/0.5


O código para ajustar o modelo quadrático e a saída são:

library(rsm)
treb.quad = rsm(y ~ SO(x1, x2, x3), data = Treb)
summary(treb.quad)
## 
## Call:
## rsm(formula = y ~ SO(x1, x2, x3), data = Treb)
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  90.00000    1.16905  76.9859 7.006e-09 ***
## x1           19.75000    0.71589  27.5880 1.171e-06 ***
## x2           19.75000    0.71589  27.5880 1.171e-06 ***
## x3          -11.50000    0.71589 -16.0639 1.703e-05 ***
## x1:x2        -6.25000    1.01242  -6.1733 0.0016247 ** 
## x1:x3         4.75000    1.01242   4.6917 0.0053768 ** 
## x2:x3         6.75000    1.01242   6.6672 0.0011461 ** 
## x1^2         -9.37500    1.05376  -8.8967 0.0002986 ***
## x2^2         -1.37500    1.05376  -1.3048 0.2487686    
## x3^2         -3.37500    1.05376  -3.2028 0.0239200 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9975, Adjusted R-squared:  0.9929 
## F-statistic: 218.9 on 9 and 5 DF,  p-value: 5.964e-06
## 
## Analysis of Variance Table
## 
## Response: y
##                 Df Sum Sq Mean Sq  F value    Pr(>F)
## FO(x1, x2, x3)   3 7299.0 2433.00 593.4146 8.448e-07
## TWI(x1, x2, x3)  3  428.8  142.92  34.8577 0.0008912
## PQ(x1, x2, x3)   3  351.5  117.16  28.5759 0.0014236
## Residuals        5   20.5    4.10                   
## Lack of fit      3   14.5    4.83   1.6111 0.4051312
## Pure error       2    6.0    3.00                   
## 
## Stationary point of response surface:
##         x1         x2         x3 
##  0.9236846 -1.7161183 -2.7698217 
## 
## Stationary point in original units:
##        A        B        C 
## 7.847369 6.419409 1.115089 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1]   1.280298  -3.551452 -11.853845
## 
## $vectors
##          [,1]       [,2]       [,3]
## x1 -0.1236692  0.5238084  0.8428112
## x2  0.8323200 -0.4077092  0.3755217
## x3  0.5403233  0.7479291 -0.3855551


A primeira tabela na saída mostra as estimativas de mínimos quadrados dos coeficientes dos fatores codificados na Equação (10.4), juntamente com seus erros padrão, valores da estatística \(t\) e \(p\)-valores. Alguns estatísticos defendem a retirada dos termos insignificantes do modelo e a reformulação, mas isso não é possível com a função rsm.

Outros estatísticos defendem a manutenção de todos os termos no modelo quadrático geral, sejam significativos ou não, uma vez que os termos do modelo individual não têm significado físico. Um compromisso é eliminar um fator do modelo se todos os termos lineares, quadráticos e de interação envolvendo esse fator forem insignificantes.

Ao eliminar um fator, o modelo de superfície de resposta é reduzido de \(p\) dimensões para \(p-1\) dimensões. Neste exemplo, todos os termos lineares e de interação são significativos, portanto nenhum fator pode ser eliminado. Se um fator fosse considerado insignificante, ele poderia ser eliminado da fórmula SO(x1, x2, x3) para a declaração do modelo quadrático e a equação da superfície de resposta seria simplificada de três dimensões para duas.

A segunda tabela mostra a tabela ANOVA do ajuste. Neste exemplo, pode-se ver que os termos lineares quadráticos e de produto cruzado no modelo são todos significativos. As somas dos quadrados dos erros são particionadas nas somas dos quadrados dos erros puros e na soma dos quadrados de bondade de ajuste. Ao fazer isso, um teste \(F\) simples pode ser feito para a adequação do modelo quadrático.

Por exemplo, as respostas nos três pontos centrais da Tabela 10.3 são 88, 91 e 91. Portanto, as somas puras dos quadrados devido à replicação são: \[ \sum_{i=1}^3 y_i^2 -\dfrac{1}{3}\Big(\sum_{i=1}^3 y_i\Big)^2 = 88^2+91^2+91^2-\dfrac{1}{3}(88+91+91)^2=6.0\cdot \]

A diferença entre as somas totais dos erros dos quadrados (20.5) e as somas dos erros puros dos quadrados é chamada de falta de ajuste das somas dos quadrados. Um teste \(F\) de falta de ajuste de somas de quadrados é um teste de adequação do modelo quadrático.

Neste caso, é claramente insignificante, indicando que as previsões feitas a partir do modelo quadrático geral para esta experiência podem ser consideradas tão precisas como a execução de experiências adicionais, desde que nenhuma variável oculta mude antes que experiências adicionais possam ser realizadas.

Quando o teste de falta de ajuste é significativo, indica que o modelo quadrático geral não é adequado para predição. Isto pode ser devido ao fato da região experimental ser tão grande que o modelo quadrático não fornece uma boa aproximação para a verdadeira função de resposta em toda a região, ou devido a outliers, ou extrema não linearidade em certos cantos da região experimental onde a função aproximada não cabe.

Essas condições podem ser detectadas fazendo gráficos de resíduos para verificar as suposições de mínimos quadrados ou ajustando o modelo usando uma técnica “robusta”, como os M-estimadores; para detalhes e discussão, consulte Lawson (1982).

Os termos de blocos também podem ser incluídos ao ajustar o modelo quadrático com a função rsm. Por exemplo, os experimentos com argamassa de cimento mostrados na Tabela 10.1 poderiam ser concluídos em dois blocos. O primeiro bloco consistiria no fatorial mais os pontos centrais mostrados nas execuções 1-11 e o segundo bloco consistiria nos pontos axiais e centrais mostrados nas execuções 12-20.

O código R para recuperar os dados e o modelo são mostrados abaixo.

library(daewr)
data(cement)
grout.quad = rsm(y ~ Block + SO(x1,x2,x3), data=cement)
summary(grout.quad)
## 
## Call:
## rsm(formula = y ~ Block + SO(x1, x2, x3), data = cement)
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)  1.1628e+02  1.0691e+00 108.7658 2.383e-15 ***
## Block2       4.4393e-01  1.0203e+00   0.4351   0.67375    
## x1           5.4068e+00  6.1057e-01   8.8553 9.746e-06 ***
## x2           9.2860e-01  6.1057e-01   1.5209   0.16262    
## x3           4.9925e+00  6.1057e-01   8.1767 1.858e-05 ***
## x1:x2        1.2500e-01  7.9775e-01   0.1567   0.87895    
## x1:x3       -1.3443e-14  7.9775e-01   0.0000   1.00000    
## x2:x3        1.2500e-01  7.9775e-01   0.1567   0.87895    
## x1^2         1.4135e+00  5.9582e-01   2.3723   0.04175 *  
## x2^2         1.3251e+00  5.9582e-01   2.2240   0.05322 .  
## x3^2         1.5019e+00  5.9582e-01   2.5207   0.03273 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9473, Adjusted R-squared:  0.8887 
## F-statistic: 16.17 on 10 and 9 DF,  p-value: 0.0001414
## 
## Analysis of Variance Table
## 
## Response: y
##                 Df Sum Sq Mean Sq F value    Pr(>F)
## Block            1   0.00   0.003  0.0006   0.98068
## FO(x1, x2, x3)   3 751.41 250.471 49.1962 6.607e-06
## TWI(x1, x2, x3)  3   0.25   0.083  0.0164   0.99693
## PQ(x1, x2, x3)   3  71.45  23.817  4.6779   0.03106
## Residuals        9  45.82   5.091                  
## Lack of fit      5  42.49   8.498 10.1972   0.02149
## Pure error       4   3.33   0.833                  
## 
## Stationary point of response surface:
##         x1         x2         x3 
## -1.9045158 -0.1825251 -1.6544845 
## 
## Stationary point in original units:
##     WatCem     BlackL        SNF 
## 0.32095484 0.14452425 0.06691031 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.525478 1.436349 1.278634
## 
## $vectors
##         [,1]       [,2]       [,3]
## x1 0.1934409  0.8924556  0.4075580
## x2 0.3466186  0.3264506 -0.8793666
## x3 0.9178432 -0.3113726  0.2461928


A saída será semelhante ao último exemplo, exceto que um termo de bloco com um grau de liberdade será incluído em ambas as tabelas. O fator de bloco será então responsável por quaisquer diferenças médias entre o primeiro e o segundo grupo de experimentos.


10.6.3 Ajustando um modelo mecanístico não linear


O modelo quadrático geral é linear nos coeficientes \(\beta\)’s e o \(\pmb{X}^\top\pmb{X}\) é de posto completo, portanto as estimativas de mínimos quadrados são obtidas resolvendo as equações normais por meio de inversão direta de matriz.

Quando o modelo é não linear como o modelo de dois compartimentos mostrado na Equação (10.7), as estimativas de mínimos quadrados são mais difíceis de obter e devem ser encontradas por técnicas numéricas iterativas. No entanto, a função R nls torna tudo isso transparente para o usuário.

Por exemplo, para ajustar o modelo de dois compartimentos (Equação 10.7) aos dados de Wagner (1967) mostrados na Tabela 10.7, use o código a seguir. A saída resultante é mostrada também a continuação. A tabela fornece as estimativas dos parâmetros de mínimos quadrados, seus erros padrão, estatísticas \(t\) e \(p\)-valores.

library(daewr)
data(Tet)
mod.nln1 = nls( Conc ~ gamma0 * (exp( -k1 * (Time - t0)) - exp( -k2 * (Time - t0))), 
                data = Tet, start = list( gamma0 = 10, k1 = 0.12, k2 = 0.5, t0 = 0.5))
summary(mod.nln1)
## 
## Formula: Conc ~ gamma0 * (exp(-k1 * (Time - t0)) - exp(-k2 * (Time - t0)))
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## gamma0  2.64964    0.36446   7.270 0.000770 ***
## k1      0.14880    0.01441  10.327 0.000146 ***
## k2      0.71575    0.12605   5.678 0.002359 ** 
## t0      0.41224    0.09495   4.342 0.007416 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04482 on 5 degrees of freedom
## 
## Number of iterations to convergence: 11 
## Achieved convergence tolerance: 3.632e-06


Neste código, os dados da Tabela 10.7 são recuperados do conjunto de dados do pacote daewr. A fórmula dada na chamada à função nls é a Equação (10.7), onde Time no conjunto de dados Tet é \(x\) na Equação (10.7). A opção start = list( gama0 = 10, k1 = 0.12, k2 = 0.5, t0 = 0.5) fornece estimativas iniciais dos parâmetros a partir dos quais o procedimento numérico iterativo começará.

A convergência do procedimento numérico pode ser sensível às estimativas iniciais, e Bates and Watts (2007) apresentam cinco estratégias para obter boas estimativas iniciais. A função nls por padrão, outras opções estão disponíveis, usa o método Gauss-Newton com derivadas numéricas para encontrar as estimativas de mínimos quadrados.


10.7 Determinação das condições operacionais ideais



10.7.1 Gráficos de contorno


Uma vez que um modelo de superfície de resposta tenha sido ajustado aos dados de um experimento, existem várias maneiras de identificar as condições ideais.

Se houver apenas dois fatores, o método mais simples é fazer um gráfico de contorno ou um gráfico tridimensional da superfície de resposta. Se houver mais de dois fatores no modelo, várias fatias bidimensionais através da região experimental podem ser feitas mantendo alguns fatores constantes e fazendo um contorno ou gráficos 3D em relação aos dois fatores restantes em cada fatia.

O pacote rsm gera automaticamente um painel de gráficos de contorno. Os comandos abaixo mostram que simplesmente adicionando a instrução contour(treb.quad, ~ x1+x2+x3) o painel de gráficos de contorno mostrado na Figura 10.15 é criado.

library(daewr)
data(Treb)
library(rsm)
treb.quad = rsm(y ~ SO(x1, x2, x3), data = Treb)
par (mar = c(4,4,1,1), mfrow=c(2,2))
contour(treb.quad, ~ x1+x2+x3 )

Figura 10.15: Gráficos de contorno da distância prevista do trabuco.

Neste painel de gráficos, o primeiro gráfico no canto superior esquerdo, é feito em uma fatia com C: peso do míssil, mantido constante em seu valor médio de 2.5 e fatores A: comprimento do braço e B: contrapeso no eixo.

As linhas representam contornos dos valores previstos do modelo quadrático. Os dois gráficos restantes mostram os fatores A e C no eixo com o fator B mantido constante em seu valor médio de 15 e os fatores B e C no eixo com A mantido constante em seu valor médio de 6.

Chamando a função gráfica persp em vez do contorno, como mostrado abaixo, é criado o painel de gráficos de superfície tridimensionais mostrado na Figura 10.16. A opção contours=list(z=“bottom”) inclui as mesmas linhas de contorno mostradas na Figura 10.15, abaixo da superfície 3D.

par (mar = c(1,1,1,1), mfrow=c(2,2))
persp(treb.quad, ~ x1+x2+x3, zlab="Distance", contours=list(z="bottom") )

Figura 10.16: Gráficos de superfície tridimensionais da distância prevista do trabuco.

Contornos individuais ou gráficos 3D também podem ser criados e os valores constantes dos fatores não mostrados no eixo também podem ser especificados conforme mostrado no código abaixo.

par (mar = c(1,1,1,1), mfrow=c(1,2))
contour(treb.quad, x1~x3, at=list(x2=1))
persp(treb.quad, x1~x3, at=list(x2=1),zlab="Distance", contours=list(z="bottom"))

Figura 10.17: Contorno e gráfico 3D da distância prevista do trabuco com contrapeso = 20 lb.

O gráfico de contorno resultante e o gráfico 3D em uma fatia com B: contrapeso = 20 lb são mostrados na Figura 10.17. Neste gráfico pode-se observar que com um contrapeso de 20 libras, a distância máxima prevista será alcançada usando um peso de míssil de cerca de 2.375 onças e um comprimento de braço de cerca de 7.5 polegadas. A distância prevista nestas condições é de cerca de 115 pés, com um erro padrão de cerca de 1.5 pés.


10.7.2 Análise canônica


Quando a superfície de resposta é um topo de colina ou um vale com um máximo ou mínimo distinto dentro da região experimental, as coordenadas exatas do fator do máximo ou mínimo podem ser determinadas definindo simultaneamente as derivadas da equação estimada em relação a cada fator igual a zero e resolver o conjunto simultâneo de equações.

Isto é útil quando há mais de dois fatores e o máximo ou mínimo seria difícil de identificar com múltiplos gráficos de contorno. A solução vetorial, de configurações de fator, para o conjunto simultâneo de equações homogêneas pode ser expressa em termos matriciais como \[\begin{equation*} \tag{10.9} \pmb{x}_0=-\widehat{\pmb{B}}^{-1}\widehat{\pmb{b}}/2 \end{equation*}\] onde \(\widehat{\pmb{B}}\) e \(\widehat{\pmb{b}}\) são as estimativas de mínimos quadrados da matriz e do vetor de coeficientes de regressão definidos na Equação (10.5). Esta solução é na verdade chamada de ponto estacionário porque pode ser um máximo, um mínimo ou um ponto de sela, conforme mostrado na Figura 10.1.

Para determinar qual é a solução, é útil expressar a Equação (10.5) da superfície de resposta em uma forma canônica com a origem traduzida para \(\pmb{x}_0\) e o eixo girado, como mostrado na Figura 10.18.

Seja \(\pmb{z} = \pmb{x}-\pmb{x}_0\), a Equação (10.5) pode ser escrita como \[\begin{equation*} \tag{10.10} \pmb{y}=\pmb{y}_0+\pmb{z}^\top \pmb{B}\pmb{z}, \end{equation*}\] onde \(\pmb{y}_0=\pmb{x}_0\pmb{b}+\pmb{x}^\top_0 \pmb{B}\pmb{x}_0\).

Isso traduz a origem para \(\pmb{x}_0\). Através da rotação ortogonal \(\pmb{z} = \pmb{M}\pmb{w}\), a parte de segunda ordem da Equação (10.10) pode ser reduzida a uma combinação linear de quadrados das variáveis rotacionadas \(w_i\) dada por \[\begin{equation*} \tag{10.11} \pmb{z}^\top\pmb{B}\pmb{z}=\pmb{w}^\top\mbox{M}^\top\mbox{BMw} = \lambda_1 w_1^2 + \lambda_2 w_2^2 + \cdots + \lambda_k w_k^2\cdot \end{equation*}\]

Se todos os coeficientes \(\lambda_i\) forem negativos, isso indica que o ponto estacionário é máximo. Se todos os coeficientes forem positivos, indica que o ponto estacionário é mínimo e se houver sinais mistos dos coeficientes, indica que o ponto estacionário é um ponto de sela. A matriz \(\pmb{M}\) é a matriz que contém os autovetores de \(\pmb{B}\) como as colunas e os coeficientes \(\lambda_i\) são os autovalores de \(\pmb{B}\).

A função rsm calcula automaticamente os autovalores e autovetores das estimativas de mínimos quadrados de \(\widehat{\pmb{B}}\). Por exemplo, para o experimento com argamassa de cimento mostrado na Tabela 10.1, o código a seguir produz a seguinte tabela de resultados.

library(daewr)
data(cement)
grout.quad = rsm(y ~ Block + SO(x1,x2,x3), data=cement)
summary(grout.quad)
## 
## Call:
## rsm(formula = y ~ Block + SO(x1, x2, x3), data = cement)
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)  1.1628e+02  1.0691e+00 108.7658 2.383e-15 ***
## Block2       4.4393e-01  1.0203e+00   0.4351   0.67375    
## x1           5.4068e+00  6.1057e-01   8.8553 9.746e-06 ***
## x2           9.2860e-01  6.1057e-01   1.5209   0.16262    
## x3           4.9925e+00  6.1057e-01   8.1767 1.858e-05 ***
## x1:x2        1.2500e-01  7.9775e-01   0.1567   0.87895    
## x1:x3       -1.3443e-14  7.9775e-01   0.0000   1.00000    
## x2:x3        1.2500e-01  7.9775e-01   0.1567   0.87895    
## x1^2         1.4135e+00  5.9582e-01   2.3723   0.04175 *  
## x2^2         1.3251e+00  5.9582e-01   2.2240   0.05322 .  
## x3^2         1.5019e+00  5.9582e-01   2.5207   0.03273 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9473, Adjusted R-squared:  0.8887 
## F-statistic: 16.17 on 10 and 9 DF,  p-value: 0.0001414
## 
## Analysis of Variance Table
## 
## Response: y
##                 Df Sum Sq Mean Sq F value    Pr(>F)
## Block            1   0.00   0.003  0.0006   0.98068
## FO(x1, x2, x3)   3 751.41 250.471 49.1962 6.607e-06
## TWI(x1, x2, x3)  3   0.25   0.083  0.0164   0.99693
## PQ(x1, x2, x3)   3  71.45  23.817  4.6779   0.03106
## Residuals        9  45.82   5.091                  
## Lack of fit      5  42.49   8.498 10.1972   0.02149
## Pure error       4   3.33   0.833                  
## 
## Stationary point of response surface:
##         x1         x2         x3 
## -1.9045158 -0.1825251 -1.6544845 
## 
## Stationary point in original units:
##     WatCem     BlackL        SNF 
## 0.32095484 0.14452425 0.06691031 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.525478 1.436349 1.278634
## 
## $vectors
##         [,1]       [,2]       [,3]
## x1 0.1934409  0.8924556  0.4075580
## x2 0.3466186  0.3264506 -0.8793666
## x3 0.9178432 -0.3113726  0.2461928


As coordenadas do ponto estacionário são fornecidas em unidades codificadas e não codificadas. A seguir, os autovalores e autovetores de \(\widehat{\pmb{B}}\) são listados.

Como todos os autovalores são positivos, o ponto estacionário é mínimo. Gráficos de contorno e gráficos de superfície 3D podem ser feitos mantendo os fatores que não estão no eixo constantes em seus valores de ponto estacionário.

Para este problema específico, não é útil conhecer as configurações dos fatores para obter uma resposta mínima. O objetivo da experimentação foi encontrar as combinações de fatores que maximizariam a trabalhabilidade da argamassa de cimento.

# slice at stationary point instead of coding value
par ( mar=c(4,4,1,1), mfrow=c(2,2))
contour(treb.quad, ~ x1+x2+x3, at = xs(treb.quad) )


Em geral, ao buscar uma resposta máxima ou mínima dentro da região experimental, o ponto estacionário só será útil se

  1. o ponto estacionário estiver dentro da região experimental e

  2. o ponto estacionário for do tipo procurado, ou seja, máximo ou mínimo.

Se a resposta estacionária estiver fora da região experimental ou se for um ponto de sela ou mínimo na busca do máximo, outro método deve ser utilizado para identificar o ótimo.


10.7.3 Análise ridge


Outro método para encontrar o ótimo dentro da região experimental é usar a análise de ridge, Hoerl (1959); Draper (1963). Este método procura encontrar o máximo ou mínimo de \[ \pmb{y} = \pmb{xb} + \pmb{x}^\top \pmb{Bx} \] sujeito à restrição de que \[ \pmb{x}^\top \pmb{x} = R^2, \] ou que as coordenadas do ótimo estejam em um raio \(R\) da origem em unidades codificadas.

A solução é obtida na ordem inversa usando multiplicadores de Lagrange. As coordenadas ótimas resultantes são as solução da equação \[\begin{equation*} \tag{10.12} \big(\pmb{B}-\mu\pmb{I}_k\big)\pmb{x} = -\pmb{b}/2\cdot \end{equation*}\]

Para resolver este sistema, escolha um valor para \(\mu\), insira-o na Equação (10.12) e resolva o vetor de configurações de fator codificado \(x\). Uma vez encontradas as configurações do fator codificado, o raio da solução é \(R =\sqrt{\sum_i x_i^2}\).

A inserção de valores de \(\mu\) maiores que o maior autovalor da matriz \(\pmb{B}\) resultará em uma solução para a resposta máxima no raio \(R\) e a inserção de valores de \(\mu\) menores que o menor autovalor da matriz \(\pmb{B}\) resultará em uma solução para a resposta mínima no raio \(R\). Tentativas e erro são necessários para encontrar a solução em um raio específico.

No entanto, a função steepest do pacote rsm faz todo o trabalho. Como exemplo, considere novamente encontrar a distância máxima que o trabuco pode lançar um objeto dentro dos alcances experimentais definidos na Tabela 10.3. Isto foi feito aproximadamente usando um gráfico de contorno na Seção 10.7.1. A análise do ponto estacionário, impressa pela função rsm para estes dados, mostrou que o ponto estacionário era um ponto de sela que estava fora da região experimental.

Portanto, o máximo dentro da região experimental deve estar no limite extremo da região experimental. Como a região experimental para um projeto de Box-Behnken, que foi usado para este experimento, é uma esfera com raio \(\sqrt{2} = 1.412\) em unidades codificadas, a chamada para a função steepest no código mostrado abaixo dá origem à lista de saída mostrado abaixo do código.

ridge = steepest(treb.quad, dist=seq(0, 1.412, by=.1), descent=FALSE)
## Path of steepest ascent from ridge analysis:
ridge
##    dist    x1    x2     x3 |     A      B      C |    yhat
## 1   0.0 0.000 0.000  0.000 | 6.000 15.000 2.5000 |  90.000
## 2   0.1 0.064 0.067 -0.037 | 6.128 15.335 2.4815 |  92.909
## 3   0.2 0.124 0.139 -0.073 | 6.248 15.695 2.4635 |  95.626
## 4   0.3 0.180 0.215 -0.105 | 6.360 16.075 2.4475 |  98.120
## 5   0.4 0.232 0.297 -0.134 | 6.464 16.485 2.4330 | 100.455
## 6   0.5 0.277 0.385 -0.158 | 6.554 16.925 2.4210 | 102.599
## 7   0.6 0.315 0.480 -0.175 | 6.630 17.400 2.4125 | 104.590
## 8   0.7 0.345 0.580 -0.185 | 6.690 17.900 2.4075 | 106.424
## 9   0.8 0.368 0.686 -0.185 | 6.736 18.430 2.4075 | 108.154
## 10  0.9 0.384 0.795 -0.177 | 6.768 18.975 2.4115 | 109.783
## 11  1.0 0.393 0.905 -0.161 | 6.786 19.525 2.4195 | 111.318
## 12  1.1 0.397 1.017 -0.137 | 6.794 20.085 2.4315 | 112.817
## 13  1.2 0.398 1.127 -0.107 | 6.796 20.635 2.4465 | 114.259
## 14  1.3 0.395 1.236 -0.073 | 6.790 21.180 2.4635 | 115.673
## 15  1.4 0.390 1.344 -0.034 | 6.780 21.720 2.4830 | 117.077


Pode-se observar que o máximo estimado em um raio \(R\) aumenta à medida que \(R\) aumenta ao longo do caminho mostrado na Figura 10.19 e que finalmente atinge um máximo de 117.077 pés no limite da região experimental, 1.412 em unidades codificadas.

Nesse ponto, os níveis dos fatores em unidades não codificadas são aproximadamente A: comprimento do braço = 6.8 polegadas; B: contrapeso = 21.8 lb; e C: peso do míssil = 2.5 onças.

Um gráfico da respostas estimadas e das coordenadas do fator em função da distância da origem, em unidades codificadas, conforme mostrado na Figura 10.20. Esta é uma forma comum de apresentar e resumir os resultados de uma análise ridge.

par (mar=c(4,4,1,1), mfrow=c(1,2))
leg.txt<-c("A","B","C")
plot(ridge$dist,ridge$yhat, type="l",xlab="radius",ylab="Max. Predicted")
grid()
plot(ridge$dist,seq(1,22,by=1.5), type="n", xlab="radius", ylab="Factors")
grid()
lines(ridge$dist,ridge$A,lty=1)
lines(ridge$dist,ridge$B,lty=2)
lines(ridge$dist,ridge$C,lty=3)
legend(1.0,19,leg.txt,lty=c(1,2,3))

Figura 10.20: Níveis de fator para atingir o valor máximo previsto.


10.7.4 Otimização não linear


A análise canônica e a análise de ridge (crista) funcionam bem para encontrar o ótimo de uma superfície de resposta se o modelo for o modelo quadrático geral e a região experimental for esférica em unidades codificadas. Se a equação da superfície de resposta não for um modelo quadrático geral ou a região experimental for irregular, o ótimo exato sempre poderá ser encontrado usando métodos numéricos mais gerais.

Como exemplo de uso de constrOptim para encontrar uma superfície de resposta ideal, considere encontrar o máximo de uma superfície de resposta mecanicista não linear. Em uma reação química, o reagente é transformado no produto com taxa constante \(\kappa_1\) e o produto se degrada decomposto em subprodutos a taxa constante \(\kappa_2\).

Da cinética de primeira ordem, a concentração do produto \(p\), no tempo \(t\) é dada pela equação \[\begin{equation*} \tag{10.13} p=-R_0 \Big(e^{-\kappa_1 t} - e^{-\kappa_2 t}\Big) \kappa_1/(\kappa_1+\kappa_2), \end{equation*}\] onde as taxas constantes são encontradas pelas equações de Arrhenius em termos de temperatura em graus Kelvin \(T\), \[ \kappa_1=a\times e^{-b\big(1/T-1/400\big)}, \quad \kappa_2 = c\times e^{-d\big(1/T-1/400\big)}\cdot \]

Se os experimentos tivessem sido conduzidos na região onde \(t = 0\) a 25 horas e \(T = 375\) a 425 graus Kelvin e os coeficientes no modelo fossem estimados em \(\widehat{R}_0 = 132.0\), \(\widehat{a} = 0.523\), \(\widehat{b} = 9847\), \(\widehat{c} = 0.20\) e \(\widehat{d} = 12327\), então o código constrOptim abaixo pode ser usado para encontrar a concentração máxima do produto.

start = c(12.5,400)
prod = function(x) {
  time = x[1]
  Temp = x[2]
  k1 = 0.523*exp(-9847*((1/Temp-1/400)))
  k2 = 0.2*exp(-12327*((1/Temp-1/400)))
  f = 132*(exp(-k1*time)-exp(-k2*time))*k1/(k1-k2) }
ui = matrix(c(1,-1,0,0,0,0,1,-1),4,2)
ci = c(0,-25,375,-425)
constrOptim(start,prod,NULL,ui,ci)
## $par
## [1]  18.15747 375.00000
## 
## $value
## [1] -82.8794
## 
## $counts
## function gradient 
##      388       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $outer.iterations
## [1] 3
## 
## $barrier.value
## [1] -0.02614004


start representa o ponto inicial no centro da região experimental. Para encontrar o máximo usando um algoritmo de minimização, foi utilizado o negativo da função codificado como prod no código acima. A função constrOptim encontra o mínimo sujeito às restrições dadas na equação matricial \[ \pmb{u}_i^\top \pmb{x}−\pmb{c}_i\geq 0\cdot \] Neste exemplo, as restrições são \(0\leq t\leq 25\) e \(375\leq T\leq 425\).

Portanto \[\begin{equation*} \tag{10.14} \pmb{u}_i=\begin{pmatrix} 1 & 0 \\ -1 & 0 \\ 0 & 1 \\ 0 & -1 \end{pmatrix}, \quad \pmb{x}=\begin{pmatrix} t \\ t \\ T \\ T \end{pmatrix}, \quad \pmb{c}_i = \begin{pmatrix} 0 \\ -25 \\ 375 \\ -425 \end{pmatrix}, \end{equation*}\] onde \(\pmb{u}_i\) e \(\pmb{c}_i\) são de definidos no código R acima.

par(mfrow=c(1,1), mar = c(4,4,1,1))
time = seq(0,25,length=52)
Temp = seq(375,425,length=52)
p = matrix(rep(0,2704),nrow=52)
for (i in 1:52)    {
  for (j in 1:52) {
    k1 = 0.523*exp(-9847*((1/Temp[j]-1/400)))
    k2 = 0.2*exp(-12327*((1/Temp[j]-1/400)))
    p[i,j] = -132*(exp(-k1*time[i])-exp(-k2*time[i]))*k1/(k1-k2)
  }
}
contour(time,Temp,p,xlab="Time, hrs",ylab="Temp, Deg K")
grid()
arrows(x0 = 20.0, y0 = 410, 
       x1 = 18.2, y1 = 375, 
       length = 0.1)

text(x = 20.0, y = 415, 
     labels = "Máximo p = 82.9")

Figura 10.21: Concentração máxima do produto.

Os resultados mostram que o mínimo do negativo da concentração do produto \(p\), dentro da região experimental é −82.88 em \(t = 18.16\) horas e \(T = 375\) graus Kelvin, o que significa que o máximo \(p\) dentro da região experimental é 82.88.

Como este problema possui apenas dois fatores (tempo e temperatura), os resultados podem ser visualizados no gráfico de contorno mostrado na Figura 10.21, que mostra o máximo no limite da região experimental.

Uma aproximação grosseira à otimização numérica pode ser feita fazendo uma pesquisa em grade. Crie uma grade fina da região experimental num onjunto de dados R e, em seguida, avalie a função da superfície de resposta em cada ponto da grade para obter valores previstos. Finalmente, classifique o conjunto de dados resultante pelos valores previstos para localizar a resposta prevista ideal dentro da região experimental.


10.7.5 Otimização de múltiplas respostas


Freqüentemente, as configurações ideais dos fatores dependerão de mais de uma resposta. Por exemplo, os regimes medicamentosos óptimos dependem da eficácia e dos efeitos secundários, a formulação óptima dos produtos alimentares depende do gosto do consumidor e dos problemas de saúde, e os processos de fabricação óptimos dependem da qualidade do produto e do custo de produção.

Uma maneira de determinar as configurações ideais dos fatores quando há mais de uma resposta é usar a otimização numérica restrita. Por exemplo, considere um problema simples com apenas duas respostas.

Os produtos de branqueamento a seco que são seguros para tecidos coloridos consistem em perborato de sódio em combinação com um ativador que aumenta a capacidade de branqueamento do perborato em água de lavagem de baixa temperatura. Experimentos de superfície de resposta foram realizados para determinar a proporção ideal de um novo ativador para perborato e a quantidade de perborato em um produto alvejante seco.

Os três fatores no estudo foram a temperatura de lavagem (Temp), a proporção de ativador para perborato (Ratio) e a quantidade de perborato (AOPPM), medida em PPM de oxigênio ativo. A faixa dos três fatores estudados foi \(70\leq \mbox{Temp} \leq 140\), \(0.5\leq \mbox{Ratio} \leq 1.5\) e \(5.0\leq \mbox{AOPPM} \leq 65.0\). A resposta foi a porcentagem de remoção de manchas de chá.

A partir dos resultados do experimento, o modelo quadrático geral, em níveis reais de fatores, relacionando a resposta aos fatores foi estimado como sendo \[ \begin{array}{rcl} \mbox{tsr} & = & -226.38 +3.375 \times \mbox{Temp}+86.5\times \mbox{Ratio}+2.646 \times \mbox{AOPPM} \\ & & -0.0128 \times \mbox{Temp}^2-17.5\times \mbox{Ratio}^2 - 0.0121\times \mbox{AOPPM}^2 -0.3857\times \mbox{Ratio}\times\mbox{Temp} \\ & & -0.0126 \times \mbox{AOPPM}\times \mbox{Temp}-0.0333 \times\mbox{AOPPM}\times\mbox{Ratio}\cdot \end{array} \]

O custo por lavagem (em centavos) foi uma função conhecida da quantidade de ativador e perborato no produto, dado pela equação \[ \mbox{cost} = 0.8313+1.27\times \mbox{Ratio}+0.37\times\mbox{Ratio}\times\mbox{AOPPM}\cdot \]

O objetivo era maximizar a remoção de manchas de chá por um custo de 10 centavos por lavagem ou menos. Este ótimo também pode ser encontrado usando a função constrOptim no R.

As restrições de contorno podem ser traduzidas para o formato \(\pmb{u}_i^\top \pmb{x}−\pmb{c}_i\geq 0\), definindo \[\begin{equation*} \tag{10.15} \pmb{u}_i=\begin{pmatrix} 1 & 0 & 0 \\ -1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & -1 \end{pmatrix}, \quad \pmb{x}=\begin{pmatrix} \mbox{Temp} \\ \mbox{Temp} \\ \mbox{Ratio} \\ \mbox{Ratio} \\ \mbox{AOPPM} \\ \mbox{AOPPM} \end{pmatrix}, \quad \pmb{c}_i = \begin{pmatrix} 70 \\ -140 \\ 0.5 \\ -1.5 \\ 5 \\ 65 \end{pmatrix}\cdot \end{equation*}\]

A restrição de custo é não linear e não pode ser incluída no \(\pmb{u}_i^\top \pmb{x}−\pmb{c}_i\geq 0\) equações de restrição linear. Em vez disso, uma função objetivo da forma: \[ −1\times \Big(\mbox{tsr}+100\times\big(\max(\mbox{cost}, 0)−10\big)\Big) \] é usada.

O código R é mostrado a continuação.

# maximize tsr for cost <=10 cents
# Note: since no R packages allow nonlinear equality constraints
# I included this constraint in the tsr function by adding the 
# 100*maximum of (cost -10) and zero 
start = c(100,.6,40)
tsrcost = function(x) {
  Temp = x[1]
  Ratio = x[2]
  AOPPM = x[3]
  tsrcost =  -( -226 + 3.375 * Temp + 86.5 * Ratio + 2.646 * AOPPM - .0128 * Temp * Temp 
                - 17.5 *Ratio *Ratio - .0121 * AOPPM * AOPPM -.3857 * Ratio * Temp - 
                  .0126 * AOPPM * Temp - .0333 * AOPPM * Ratio) + 
    100 * max(.8313 + 1.27 * Ratio + .37 * Ratio * AOPPM - 10, 0) }
ui = matrix(c(1, -1, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 1, -1) ,6 ,3 )
ci = c(70, -140, .5, -1.5, 5, -65 ) 
constrOptim(start, tsrcost, NULL, ui, ci)
## $par
## [1] 102.9643645   0.5460873  41.9454227
## 
## $value
## [1] -40.6523
## 
## $counts
## function gradient 
##      440       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $outer.iterations
## [1] 3
## 
## $barrier.value
## [1] -0.04545005
# note this minimizes tsrcost, subject to ui%*%x-ci>=0


A execução deste código leva \(\big(\max(\mbox{cost}, 0)-10\big)\) a zero e mostra que a porcentagem máxima de remoção de manchas de chá que pode ser esperada para o ativador que está sendo testado é de 40.65% em \(\mbox{Temp = 103}\) graus, \(\mbox{Ratio = 0.55}\) e \(\mbox{AOPPM = 41.9}\).

Quando há múltiplas respostas \(\kappa\), outro método de otimização que às vezes é útil é combinar as respostas em uma função de desejabilidade proposta por Derringer and Suich (1980). A ideia é converter cada resposta em uma função de desejabilidade individual \(0\leq d_i(y_i)\leq 1\), onde esta função é 1 se a resposta \(y_i\) obtiver seu valor mais desejável e zero quando a resposta estiver em uma faixa inaceitável.

A seguir, a função de desejabilidade geral é definida como \[\begin{equation*} \tag{10.16} \big(d_1\times d_2\times \cdots\times d_\kappa\big)^{(1/\kappa)} \end{equation*}\] e esta função pode ser maximizada usando otimização numérica ou pesquisa em grade.

Se for desejável obter um valor alvo para uma resposta \(\widehat{y}_i\), por exemplo, a função de desejabilidade individual pode ser definida como \[\begin{equation*} \tag{10.17} d_i = \left\{ \begin{array}{ccc} 0 & \mbox{caso} & \widehat{y}_i<L \\ \displaystyle \left(\dfrac{\widehat{y}_i-L}{T-L}\right)^r & \mbox{caso} & L\leq \widehat{y}_i\leq T \\ \displaystyle \left(\dfrac{\widehat{y}_i-U}{T-U}\right)^s & \mbox{caso} & T\leq \widehat{y}_i\leq U \\ 0 & \mbox{caso} & \widehat{y}_i>U \end{array} \right. \end{equation*}\] onde \(L\leq T\leq U\), \(T\) é o valor alvo para a resposta e as regiões onde \(\widehat{y}_i < L\) ou \(\widehat{y}_i > U\) são indesejáveis.

As potências \(r\) e \(s\) controlam o quão crítico o pesquisador sente que é estar perto do alvo. Se for desejável que uma das respostas seja maximizada, modifique a Equação (10.15) estabelecendo \(T = U\) e definindo \(d_i = 1\) para \(\widehat{y}_i\geq U\) e \(d_i = (\widehat{y}_i − L)/(U-L)^s\) se \(L\leq y_i \leq U\). Se for desejável minimizar uma resposta, modifique a Equação (10.15) definindo \(T = L\) e definindo \(d_i = 1\) para \(\widehat{y}_i\leq L\) e \(d_i = (U-\widehat{y}_i)/(U-L)^s\) se \(L\leq \widehat{y}_i\leq U\).

A pacote R desirability (Kuhn (2013)) possui as funções dTarget, dMax e dMin para as funções de desejabilidade alvo, máximo e mínimo definidas na Equação (10.17) e a função dOverall para a função de desejabilidade geral dada na Equação (10.16). Eles podem ser usados para maximizar a Equação (10.16).

Por exemplo, no processo de fabricação de um suporte de catalisador, foi conduzido um experimento de superfície de resposta de três níveis e equações quadráticas gerais foram ajustadas para três das características do produto resultantes. Estas equações são fornecidas abaixo, onde \(x_1\), \(x_2\) e \(x_3\) são valores codificados dos fatores: tempo de mistura, tempo de filtragem e densidade de empacotamento antes da calcinação. \[ \begin{array}{rcl} \mbox{Área da Superfície} & = & 125.4106 - 8.1233 x_1 +17.0266 x_2 +0.4277 x_3 +33.88054 x_2^2 \\ & & +14.81976 x_2^2 +13.070 x_3^2+2.4184 x_1 x_2 -8.4376 x_1 x_3 +9.0134 x_2 x_3, \\ & & \\ \mbox{Volume do poro} & = & 0.661354 - 0.1963 x_1 - 0.02016 x_2 - 0.00291 x_3 + 0.15126 x^2_1 \\ & & + 0.118423 x_2^2 + 0.0679 x^2_3 + 0.02399 x_1 x_2 + 0.010327 x_1 x_3 - 0.0374 x_2 x_3, \\ & & \\ \mbox{Diâmetro do poro} & = & 39.35608 + 3.19547 x_1 + 0.2173 x_2 - 1.4698 x_3 + 0.4141 x^2_1 \\ & & - 2.39408 x_2^2 - 2.36399 x^2_3 + 0.5887 x_1 x_2 - 0.62136 x_1 x_3 - 1.53234 x_2 x_3\cdot \end{array} \]

Uma vez desenvolvidas essas equações, elas foram utilizadas para identificar as condições de processamento que produzirão um suporte catalítico com características adequadas para aplicações específicas. Por exemplo, para produzir um suporte de catalisador com Diâmetro de Poro \(\approx\) 40, Área de Superfície > 100 e Volume de Poro >0.6, as funções de desejabilidade foram de definidas conforme mostrado no código R a seguir.

O máximo da função de desejabilidade foi usado para área superficial com \(L = 100\) e \(U = 217\). O máximo da função de desejabilidade foi usado para o volume de poro com \(L = 0.6\) e \(U = 1.3\), e o alvo é melhor com a função de desejabilidade usada para diâmetro dos poros com \(L = 38\), \(T = 40\) e \(U = 42\).

x = c(0, 0, 0)
saPred = function(x) {
  125.4106 -8.1233 * x[1] +17.0266 * x[2] +.4277 * x[3] +2.4184 * x[1] * x[2] 
  -8.4376 * x[1] * x[3] +9.0134 * x[2] * x[3] + 33.88054 * x[1]^2 +14.81976 * x[2]^2 
  +13.07001 * x[3]^2}
pvPred = function(x) { .661354 -.1963 * x[1] -.02016 * x[2] -.00291 * x[3] 
  +.02399 * x[1] * x[2] +.010327 * x[1] * x[3] -.0374 * x[2] * x[3] +.15126 *x[1]^2 
  + .118423 * x[2]^2 +.0679*x[3]^2}
dpPred = function(x) { 39.35608 + 3.19547 * x[1] + .21729 * x[2] -1.46979 * x[3] + 
  .58873 * x[1] * x[2] -.62136 * x[1] * x[3] -1.53234 * x[2] * x[3] +.41413 * x[1]^2 
  -2.39408 * x[2]^2 -2.36399 * x[3]^2}
library(desirability)
saD = dMax(100, 217)
pvD = dMax(0.6, 1.3)
dpD = dTarget(38, 40, 42)
overallD = dOverall(saD, pvD, dpD)
# Code on web page referred to on p. 426
rsmOpt = function(x, dObject, space = "square")
{
  sa = saPred(x)
  pv = pvPred(x)
  dp = dpPred(x)
  
  out = predict(dObject, data.frame(sa = sa, pv = pv, dp = dp))
  
  if(space == "circular")
  {
    if(sqrt(sum(x^2)) > 1.0) out <- 0
  } else if(space == "square") if(any(abs(x) > 1.0)) out <- 0
  out
}
searchGrid = expand.grid(Mixtime = seq(-1.0, 1.0, length = 5),
                         filTtime = seq(-1.0, 1.0, length = 5),
                         PackMth = seq(-1.0, 1.0, length = 5))
for(i in 1:dim(searchGrid)[1])
{
  tmp = optim(as.vector(searchGrid[i,]),
               rsmOpt,
               dObject = overallD,
               space = "square",
               control = list(fnscale = -1))
  if(i == 1)
  {
    best <- tmp
  } else {
    if(tmp$value > best$value) best <- tmp
  }
}
predOutcomes = c(saPred(c(-0.2924993,-1,-1)), pvPred(c(-0.2924993,-1,-1)), dpPred(c(-0.2924993,-1,-1)))
print(predOutcomes)
## [1] 13.070010  0.186323 -4.758070


Depois de maximizar o overallD, os resultados indicam que a área de superfície = 111.04, o volume de poro = 0.75 e o diâmetro de poro = 39.85 podem ser alcançados quando o tempo de mistura está próximo da faixa intermediária, ou seja, \(x_1 = -0.2925\), o tempo de filtragem está no limite superior da faixa testada, ou seja, \(x_2 = 1.0\) e a densidade de empacotamento está na extremidade inferior da faixa testada, ou seja, \(x_3 = −1.0\).


10.8 Projetos em blocos de superfície de resposta


Quando as unidades experimentais não são homogêneas é sempre vantajoso agrupá-las em blocos mais homogêneos e utilizar um delineamento experimental em blocos. Anteriormente foi mencionado que os experimentos compostos centrais de três fatores poderiam ser executados em dois blocos; a porção fatorial mais três pontos centrais em um bloco e a porção axial e os demais pontos centrais no outro bloco.

Na verdade, este é um exemplo de projeto de blocos incompletos, como os descritos no Capítulo 7, uma vez que nem todas as combinações de tratamento estão representadas em cada bloco. É possível criar um projeto de bloco incompleto a partir de qualquer projeto de superfície de resposta para que os coeficientes possam ser estimados com mais precisão.

Projetos de superfície de resposta padrão, como os projetos compostos centrais (CCDs) e projetos Box-Behnken (BBDs) foram planejados em blcos, de forma que os blocos sejam ortogonais aos níveis de fator codificados, quadrados de níveis de fator codificados e interações entre níveis de fator codificados.

Desta forma, as estimativas de mínimos quadrados dos parâmetros do modelo quadrático geral não estão correlacionadas com os efeitos de bloco e estes planejamentos padrão são 100% \(D_s\) eficientes, conforme descrito na Seção 7.8.1.

O exemplo descrito anteriormente, onde o experimento composto central de três fatores foi planejado em dois blocos com tamanhos de bloco 11 e 9, não é 100% \(D_s\) eficiente. Para atingir 100% de eficiência, algumas restrições adicionais devem ser feitas em termos do número de blocos, dos tamanhos dos blocos e do raio axial.

Isto afetará a distribuição da variância de um valor previsto na região de projeto e um compromisso deve ser feito. A Tabela 10.8 mostra o número de blocos e tamanhos de bloco para CCDs e BBDs planejados em blocos ortogonais com \(k = 2\) a 5 fatores.

\[ \begin{array}{c|cc|cc}\hline \mbox{Número de} & \mbox{CCD} & & \mbox{BBD} & \\ \mbox{fatores} & \mbox{No. de blocos} & \mbox{Tamanho dos blocos} & \mbox{No. de blocos} & \mbox{Tamanho dos blocos} \\\hline 2 & 2 & 7, 7 & & & \\ 3 & 3 & 6, 6, 8 & & & \\ 4 & 3 & 10, 10, 10 & 3 & 9, 9, 9 \\ 5 & 2 & 22, 11 & 2 & 23, 23 \\\hline \end{array} \] Tabela 10.8: Número de blocos e tamanhos de blocos para projetos CCD e BBD planejados em blocos ortogonais.

Esses designs podem ser criados com as funções bbd e ccd no pacote R rsm. Por exemplo, o código para criar os experimentos Box-Behnken planejados em blocos ortogonais com quatro e cinco fatores listados na Tabela 10.8 é mostrado abaixo.

O número de pontos centrais \(n_0\) deve ser escolhido de modo que os tamanhos dos blocos correspondam aos listados na Tabela 10.8.

library(rsm)
bbd(4,n0=1,randomize=FALSE)
##    run.order std.order Block x1.as.is x2.as.is x3.as.is x4.as.is
## 1          1         1     1       -1       -1        0        0
## 2          2         2     1        1       -1        0        0
## 3          3         3     1       -1        1        0        0
## 4          4         4     1        1        1        0        0
## 5          5         5     1        0        0       -1       -1
## 6          6         6     1        0        0        1       -1
## 7          7         7     1        0        0       -1        1
## 8          8         8     1        0        0        1        1
## 9          9         9     1        0        0        0        0
## 10        10        10     2       -1        0        0       -1
## 11        11        11     2        1        0        0       -1
## 12        12        12     2       -1        0        0        1
## 13        13        13     2        1        0        0        1
## 14        14        14     2        0       -1       -1        0
## 15        15        15     2        0        1       -1        0
## 16        16        16     2        0       -1        1        0
## 17        17        17     2        0        1        1        0
## 18        18        18     2        0        0        0        0
## 19        19        19     3       -1        0       -1        0
## 20        20        20     3        1        0       -1        0
## 21        21        21     3       -1        0        1        0
## 22        22        22     3        1        0        1        0
## 23        23        23     3        0       -1        0       -1
## 24        24        24     3        0        1        0       -1
## 25        25        25     3        0       -1        0        1
## 26        26        26     3        0        1        0        1
## 27        27        27     3        0        0        0        0
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is
## x4 ~ x4.as.is
bbd(5,n0=3,randomize=FALSE)
##    run.order std.order Block x1.as.is x2.as.is x3.as.is x4.as.is x5.as.is
## 1          1         1     1       -1       -1        0        0        0
## 2          2         2     1        1       -1        0        0        0
## 3          3         3     1       -1        1        0        0        0
## 4          4         4     1        1        1        0        0        0
## 5          5         5     1       -1        0       -1        0        0
## 6          6         6     1        1        0       -1        0        0
## 7          7         7     1       -1        0        1        0        0
## 8          8         8     1        1        0        1        0        0
## 9          9         9     1        0        0       -1       -1        0
## 10        10        10     1        0        0        1       -1        0
## 11        11        11     1        0        0       -1        1        0
## 12        12        12     1        0        0        1        1        0
## 13        13        13     1        0        0        0       -1       -1
## 14        14        14     1        0        0        0        1       -1
## 15        15        15     1        0        0        0       -1        1
## 16        16        16     1        0        0        0        1        1
## 17        17        17     1        0       -1        0        0       -1
## 18        18        18     1        0        1        0        0       -1
## 19        19        19     1        0       -1        0        0        1
## 20        20        20     1        0        1        0        0        1
## 21        21        21     1        0        0        0        0        0
## 22        22        22     1        0        0        0        0        0
## 23        23        23     1        0        0        0        0        0
## 24        24        24     2       -1        0        0       -1        0
## 25        25        25     2        1        0        0       -1        0
## 26        26        26     2       -1        0        0        1        0
## 27        27        27     2        1        0        0        1        0
## 28        28        28     2       -1        0        0        0       -1
## 29        29        29     2        1        0        0        0       -1
## 30        30        30     2       -1        0        0        0        1
## 31        31        31     2        1        0        0        0        1
## 32        32        32     2        0       -1       -1        0        0
## 33        33        33     2        0        1       -1        0        0
## 34        34        34     2        0       -1        1        0        0
## 35        35        35     2        0        1        1        0        0
## 36        36        36     2        0       -1        0       -1        0
## 37        37        37     2        0        1        0       -1        0
## 38        38        38     2        0       -1        0        1        0
## 39        39        39     2        0        1        0        1        0
## 40        40        40     2        0        0       -1        0       -1
## 41        41        41     2        0        0        1        0       -1
## 42        42        42     2        0        0       -1        0        1
## 43        43        43     2        0        0        1        0        1
## 44        44        44     2        0        0        0        0        0
## 45        45        45     2        0        0        0        0        0
## 46        46        46     2        0        0        0        0        0
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is
## x4 ~ x4.as.is
## x5 ~ x5.as.is


Para criar os experimentos compostos centrais planejados em blocos ortogonais de 2 a 5 fatores, o código a seguir é usado.

library(rsm)
ccd(2, n0 = c(3, 3), alpha = "orthogonal", randomize = FALSE)
##    run.order std.order  x1.as.is  x2.as.is Block
## 1          1         1 -1.000000 -1.000000     1
## 2          2         2  1.000000 -1.000000     1
## 3          3         3 -1.000000  1.000000     1
## 4          4         4  1.000000  1.000000     1
## 5          5         5  0.000000  0.000000     1
## 6          6         6  0.000000  0.000000     1
## 7          7         7  0.000000  0.000000     1
## 8          1         1 -1.414214  0.000000     2
## 9          2         2  1.414214  0.000000     2
## 10         3         3  0.000000 -1.414214     2
## 11         4         4  0.000000  1.414214     2
## 12         5         5  0.000000  0.000000     2
## 13         6         6  0.000000  0.000000     2
## 14         7         7  0.000000  0.000000     2
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
ccd(3, n0 = 2, alpha = "orthogonal", randomize = FALSE, blocks = Block ~ (x1*x2*x3))
##     run.order std.order  x1.as.is  x2.as.is  x3.as.is Block
## 1           1         1 -1.000000 -1.000000 -1.000000     1
## 4           2         2  1.000000  1.000000 -1.000000     1
## 6           3         3  1.000000 -1.000000  1.000000     1
## 7           4         4 -1.000000  1.000000  1.000000     1
## 11          5         5  0.000000  0.000000  0.000000     1
## 2           6         6  0.000000  0.000000  0.000000     1
## 12          1         1  1.000000 -1.000000 -1.000000     2
## 41          2         2 -1.000000  1.000000 -1.000000     2
## 61          3         3 -1.000000 -1.000000  1.000000     2
## 71          4         4  1.000000  1.000000  1.000000     2
## 111         5         5  0.000000  0.000000  0.000000     2
## 21          6         6  0.000000  0.000000  0.000000     2
## 13          1         1 -1.632993  0.000000  0.000000     3
## 22          2         2  1.632993  0.000000  0.000000     3
## 3           3         3  0.000000 -1.632993  0.000000     3
## 42          4         4  0.000000  1.632993  0.000000     3
## 5           5         5  0.000000  0.000000 -1.632993     3
## 62          6         6  0.000000  0.000000  1.632993     3
## 72          7         7  0.000000  0.000000  0.000000     3
## 8           8         8  0.000000  0.000000  0.000000     3
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is
ccd(4, n0 = 2, alpha = "orthogonal", randomize = FALSE, blocks = Block ~ (x1*x2*x3*x4))
##     run.order std.order x1.as.is x2.as.is x3.as.is x4.as.is Block
## 2           1         1        1       -1       -1       -1     1
## 3           2         2       -1        1       -1       -1     1
## 5           3         3       -1       -1        1       -1     1
## 8           4         4        1        1        1       -1     1
## 9           5         5       -1       -1       -1        1     1
## 12          6         6        1        1       -1        1     1
## 14          7         7        1       -1        1        1     1
## 15          8         8       -1        1        1        1     1
## 1           9         9        0        0        0        0     1
## 21         10        10        0        0        0        0     1
## 22          1         1       -1       -1       -1       -1     2
## 31          2         2        1        1       -1       -1     2
## 51          3         3        1       -1        1       -1     2
## 81          4         4       -1        1        1       -1     2
## 91          5         5        1       -1       -1        1     2
## 121         6         6       -1        1       -1        1     2
## 141         7         7       -1       -1        1        1     2
## 151         8         8        1        1        1        1     2
## 11          9         9        0        0        0        0     2
## 211        10        10        0        0        0        0     2
## 13          1         1       -2        0        0        0     3
## 23          2         2        2        0        0        0     3
## 32          3         3        0       -2        0        0     3
## 4           4         4        0        2        0        0     3
## 52          5         5        0        0       -2        0     3
## 6           6         6        0        0        2        0     3
## 7           7         7        0        0        0       -2     3
## 82          8         8        0        0        0        2     3
## 92          9         9        0        0        0        0     3
## 10         10        10        0        0        0        0     3
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is
## x4 ~ x4.as.is
ccd(4, generators = (x5 ~ x1*x2*x3*x4), n0 = c(6,1), alpha = "orthogonal", randomize = FALSE)
##    run.order std.order x1.as.is x2.as.is x3.as.is x4.as.is x5.as.is Block
## 1          1         1       -1       -1       -1       -1        1     1
## 2          2         2        1       -1       -1       -1       -1     1
## 3          3         3       -1        1       -1       -1       -1     1
## 4          4         4        1        1       -1       -1        1     1
## 5          5         5       -1       -1        1       -1       -1     1
## 6          6         6        1       -1        1       -1        1     1
## 7          7         7       -1        1        1       -1        1     1
## 8          8         8        1        1        1       -1       -1     1
## 9          9         9       -1       -1       -1        1       -1     1
## 10        10        10        1       -1       -1        1        1     1
## 11        11        11       -1        1       -1        1        1     1
## 12        12        12        1        1       -1        1       -1     1
## 13        13        13       -1       -1        1        1        1     1
## 14        14        14        1       -1        1        1       -1     1
## 15        15        15       -1        1        1        1       -1     1
## 16        16        16        1        1        1        1        1     1
## 17        17        17        0        0        0        0        0     1
## 18        18        18        0        0        0        0        0     1
## 19        19        19        0        0        0        0        0     1
## 20        20        20        0        0        0        0        0     1
## 21        21        21        0        0        0        0        0     1
## 22        22        22        0        0        0        0        0     1
## 23         1         1       -2        0        0        0        0     2
## 24         2         2        2        0        0        0        0     2
## 25         3         3        0       -2        0        0        0     2
## 26         4         4        0        2        0        0        0     2
## 27         5         5        0        0       -2        0        0     2
## 28         6         6        0        0        2        0        0     2
## 29         7         7        0        0        0       -2        0     2
## 30         8         8        0        0        0        2        0     2
## 31         9         9        0        0        0        0       -2     2
## 32        10        10        0        0        0        0        2     2
## 33        11        11        0        0        0        0        0     2
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is
## x4 ~ x4.as.is
## x5 ~ x5.as.is


Para experimentos compostos centrais (CCD), não apenas os tamanhos dos blocos devem corresponder à Tabela 10.8 para planos em blocos ortogonais, mas também a opção alpha = “orthogonal” deve ser usada.

Para os experimentos de três e quatro fatores, a parcela fatorial deve ser dividida em dois blocos utilizando a opção blocks = Block - (x1*x2*x3) e blocks = Block ~ (x1*x2*x3*x4). Para o projeto de 5 fatores, a porção fatorial é na verdade uma meia fração com os geradores generators = (x5 ~ x1*x2*x3*x4).

O número de blocos e tamanhos de blocos disponíveis para projetos de CCD e BBD planejados em blocos ortogonais pode ser muito restritivo para algumas aplicações. Por exemplo, Gilmour and Trinca (2000) apresentaram os dados de um experimento para determinar como as características da massa folhada dependiam dos parâmetros do processo de mistura.

No experimento foram variados três fatores, a vazão de alimentação (\(FR\)), o teor de umidade inicial (\(MC\)) e a velocidade da rosca (\(SS\)). Uma das respostas (\(y\)) mediu a refletância da luz da massa resultante em faixas específicas do espectro para ver como a cor da massa foi afetada pelos parâmetros de mistura. Apenas quatro execuções ou experimentos poderiam ser realizados em um dia e os experimentadores anteciparam a variação diária.

Nenhum dos projetos da Tabela 10.8 possui blocos de tamanho quatro. Assim, os experimentadores executaram o projeto modificado mostrado na Tabela 10.9, com os valores dos fatores codificados \[ x_1 = (FR −37.5)/7.5, \qquad x_2 = (MC - 21)/3 \quad \mbox{e} \quad x_3 = (SS −350)/50\cdot \]

\[ \begin{array}{ccccc}\hline \mbox{Bloco} & x_1 & x_2 & x_3 & y \\\hline 1 & -1 & -1 & -1 & 12.92 \\ 1 & -1 & 1 & 1 & 13.91 \\ 1 & 1 & -1 & 1 & 11.66 \\ 1 & 1 & 1 & -1 & 14.48 \\ 2 & -1 & -1 & 1 & 10.76 \\ 2 & -1 & 1 & -1 & 14.41 \\ 2 & 1 & -1 & -1 & 12.27 \\ 2 & 1 & 1 & 1 & 12.13 \\ 3 & -1 & 1 & -1 & 14.22 \\ 3 & 0 & -1 & 0 & 12.35 \\ 3 & 1 & 0 & 0 & 13.50 \\ 3 & 0 & 0 & 1 & 12.54 \\ 4 & 1 & -1 & 1 & 10.55 \\ 4 & -1 & 0 & 0 & 13.33 \\ 4 & 0 & 1 & 0 & 13.84 \\ 4 & 0 & 0 & -1 & 14.19 \\ 5 & -1 & -1 & -1 & 11.46 \\ 5 & 1 & 1 & 1 & 11.32 \\ 5 & 0 & 0 & 0 & 11.93 \\ 5 & 0 & 0 & 0 & 11.63 \\ 6 & -1 & -1 & 1 & 12.20 \\ 6 & 1 & 1 & -1 & 14.78 \\ 6 & 0 & 0 & 0 & 14.94 \\ 6 & 0 & 0 & 0 & 14.61 \\ 7 & -1 & 1 & 1 & 12.17 \\ 7 & 1 & -1 & -1 & 11.28 \\ 7 & 0 & 0 & 0 & 11.85 \\ 7 & 0 & 0 & 0 & 11.64 \\\hline \end{array} \] Tabela 10.9: Projeto e resposta do experimento de massa de pastelaria.

Este é um exemplo de projeto de cubo centrado na face que é um caso especial de projeto composto central com os pontos axiais, em unidades codificadas, puxados para a face do cubo a \(\pm\) 1.

A porção fatorial \(2^3\) deste planejamento foi replicada e havia seis pontos centrais. As 28 execuções deste desenho foram pensadas em sete blocos (dias) de quatro execuções de forma que os efeitos principais sejam ortogonais aos blocos. No entanto, como apontou Goos (2002), os termos quadráticos e de interação do modelo quadrático geral não são ortogonais aos blocos e uma melhor maneira planejar os blocos de 28 experimentos pode ser alcançada usando um planejamento \(D_s\)-ótimo. Isso pode ser feito usando o pacote AlgDesign.

Primeiro crie o plano do cubo centrado na face usando a função gen.factorial e o código R mostrado abaixo. Isso cria todos os 28 pontos de projeto mostrados na Tabela 10.9 e os armazena no conjunto de dados cand com os fatores rotulados como X1, X2 e X3.

library(AlgDesign)
fact = gen.factorial(levels = 2,nVars = 3)
fact = rbind(fact, fact)
center = data.frame( matrix( rep (c(0, 0, 0), 6),ncol = 3))
star = data.frame( rbind( diag(3), -diag(3) ))
cand = rbind(fact, center, star)
bdesign = optBlock (~ quad(.), cand, blocksizes = c(4, 4, 4, 4, 4, 4, 4), 
                     criterion = "Dpc", nRepeats = 1000)
bdesign
## $D
## [1] 0.3028924
## 
## $Dpc
## [1] 0.3047534
## 
## $Blocks
## $Blocks$B1
##    X1 X2 X3
## 10  1 -1 -1
## 11 -1  1 -1
## 12  1  1 -1
## 17  0  0  0
## 
## $Blocks$B2
##    X1 X2 X3
## 2   1 -1 -1
## 5  -1 -1  1
## 8   1  1  1
## 24  0  1  0
## 
## $Blocks$B3
##    X1 X2 X3
## 13 -1 -1  1
## 21  0  0  0
## 25  0  0  1
## 26 -1  0  0
## 
## $Blocks$B4
##    X1 X2 X3
## 7  -1  1  1
## 19  0  0  0
## 23  1  0  0
## 27  0 -1  0
## 
## $Blocks$B5
##    X1 X2 X3
## 9  -1 -1 -1
## 14  1 -1  1
## 20  0  0  0
## 28  0  0 -1
## 
## $Blocks$B6
##    X1 X2 X3
## 1  -1 -1 -1
## 3  -1  1 -1
## 15 -1  1  1
## 18  0  0  0
## 
## $Blocks$B7
##    X1 X2 X3
## 4   1  1 -1
## 6   1 -1  1
## 16  1  1  1
## 22  0  0  0
## 
## 
## $design
##    X1 X2 X3
## 10  1 -1 -1
## 11 -1  1 -1
## 12  1  1 -1
## 17  0  0  0
## 2   1 -1 -1
## 5  -1 -1  1
## 8   1  1  1
## 24  0  1  0
## 13 -1 -1  1
## 21  0  0  0
## 25  0  0  1
## 26 -1  0  0
## 7  -1  1  1
## 19  0  0  0
## 23  1  0  0
## 27  0 -1  0
## 9  -1 -1 -1
## 14  1 -1  1
## 20  0  0  0
## 28  0  0 -1
## 1  -1 -1 -1
## 3  -1  1 -1
## 15 -1  1  1
## 18  0  0  0
## 4   1  1 -1
## 6   1 -1  1
## 16  1  1  1
## 22  0  0  0
## 
## $rows
##  [1] 10 11 12 17  2  5  8 24 13 21 25 26  7 19 23 27  9 14 20 28  1  3 15 18  4
## [26]  6 16 22


Em seguida, a função optBlock é usada para agrupar as execuções em sete blocos de quatro execuções. Isso mantém todas as 28 execuções no conjunto de dados e as organiza em sete blocos de quatro de forma a maximizar a eficiência do \(D_s\).

Também pode ser criado um projeto de superfície de resposta em blocos que seja um subconjunto das execuções possíveis em um arquivo de candidatos. Por exemplo, para planejar em 20 blocos os experimentos no arquivo candidato QSAR discutido na Seção 10.5 para que possam ser executados em cinco laboratórios diferentes, que podem ter diferentes medições de bioatividade, modifique o código da Seção 10.5 para aquele mostrado abaixo.

library(daewr)
data(qsar)
library(AlgDesign)
desgn1 = optBlock( ~ quad(.), qsar, blocksizes = c(4, 4, 4, 4, 4), 
                   criterion = "Dpc", nRepeats = 1000)


Para analisar um projeto em blocos com mais de dois blocos usando a função rsm, certifique-se de que o indicador Block seja do fator de classe conforme mostrado no código abaixo.

library(daewr)
data(pastry)
class(pastry$Block)
## [1] "factor"
library(rsm)
blkrsm = rsm(y ~ Block + SO(x1, x2, x3), data = pastry) 
summary(blkrsm)
## 
## Call:
## rsm(formula = y ~ Block + SO(x1, x2, x3), data = pastry)
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 13.952045   0.224989 62.0122 < 2.2e-16 ***
## Block2      -0.850000   0.220043 -3.8629  0.002257 ** 
## Block3      -0.432828   0.237417 -1.8231  0.093286 .  
## Block4      -0.607828   0.237417 -2.5602  0.024993 *  
## Block5      -1.976069   0.246980 -8.0009 3.756e-06 ***
## Block6       0.688931   0.246980  2.7894  0.016362 *  
## Block7      -2.076069   0.246980 -8.4058 2.258e-06 ***
## x1          -0.189444   0.073348 -2.5828  0.023973 *  
## x2           0.878333   0.073348 11.9749 4.950e-08 ***
## x3          -0.709444   0.073348 -9.6723 5.127e-07 ***
## x1:x2       -0.189907   0.088153 -2.1543  0.052241 .  
## x1:x3       -0.060093   0.088153 -0.6817  0.508382    
## x2:x3        0.177593   0.088153  2.0146  0.066918 .  
## x1^2        -0.113182   0.187654 -0.6031  0.557642    
## x2^2        -0.433182   0.187654 -2.3084  0.039590 *  
## x3^2        -0.163182   0.187654 -0.8696  0.401583    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9753, Adjusted R-squared:  0.9445 
## F-statistic: 31.63 on 15 and 12 DF,  p-value: 2.491e-07
## 
## Analysis of Variance Table
## 
## Response: y
##                 Df  Sum Sq Mean Sq F value    Pr(>F)
## Block            6 19.5309  3.2552 33.6144 7.957e-07
## FO(x1, x2, x3)   3 23.5921  7.8640 81.2079 3.075e-08
## TWI(x1, x2, x3)  3  0.8557  0.2852  2.9455  0.075964
## PQ(x1, x2, x3)   3  1.9645  0.6548  6.7623  0.006378
## Residuals       12  1.1621  0.0968                  
## Lack of fit      5  0.6403  0.1281  1.7183  0.248424
## Pure error       7  0.5217  0.0745                  
## 
## Stationary point of response surface:
##        x1        x2        x3 
## -1.333065  1.025086 -1.370525 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.0569490 -0.1738053 -0.4787912
## 
## $vectors
##          [,1]       [,2]       [,3]
## x1  0.8067515 0.54639939 -0.2249439
## x2 -0.3207343 0.08520655 -0.9433289
## x3 -0.4962676 0.83317918  0.2439895

10.9 Projetos Split-Plot de superfície de resposta


Quando os projetos de superfície de resposta contêm fatores de diferentes etapas do processo ou alguns fatores são difíceis de variar, pode ser caro ou inconveniente randomizar completamente as execuções em um projeto de superfície de resposta padrão.

Por exemplo, considere um planejamento de superfície de resposta simples para identificar a temperatura e o tempo de cozimento que resultariam no bolo úmido ideal. O projeto e os dados são mostrados na Tabela 10.10, onde \[ x_1 = (\mbox{temperatura de cozimento} - 350)/25, \]

\[ x_2 = (\mbox{tempo de cozimento} - 31.5)/4, \] e \(y\) é a leitura de um testador de umidade.

\[ \begin{array}{cccc}\hline \mbox{Forno funcionando} & x_1 & x_2 & x_3 \\\hline 1 & -1 & -1 & 2.7 \\ 1 & -1 & 1 & 2.5 \\ 1 & -1 & 0 & 2.7 \\ 2 & 1 & -1 & 2.9 \\ 2 & 1 & 1 & 1.3 \\ 2 & 1 & 0 & 2.2 \\ 3 & 0 & -1 & 3.7 \\ 3 & 0 & 1 & 2.9 \\ 4 & 0 & 0 & 2.9 \\ 4 & 0 & 0 & 2.8 \\ 4 & 0 & 0 & 2.9 \\\hline \end{array} \] Tabela 10.10: Dados para o experimento de cozimento de bolo.

O plano é um experimento de cubo centrado na face com três níveis para cada fator e três pontos centrais. Os experimentos não foram realizados em ordem aleatória. Depois que a temperatura do forno foi ajustada e o forno pré-aquecido, até quatro bolos poderiam ser assados de uma vez.

A realização dos experimentos em uma ordem completamente aleatória exigiria 11 ciclos de forno, mas ao assar até três bolos ao mesmo tempo, como mostrado na tabela, foram necessários apenas quatro ciclos de forno. As corridas no forno foram realizadas em ordem aleatória com os dois ou três bolos de cada corrida colocados em posições aleatórias no forno.

Os bolos em cada passagem foram retirados do forno em tempos de cozedura diferentes, excepto para a passagem 4 no forno, quando todos os bolos foram cozidos durante 31.5 minutos. Ao executar os experimentos dessa maneira, resulta um projeto do tipo gráfico dividido de superfície de resposta (RSSP), onde \(x_1\), a temperatura de cozimento, é o fator de todo o gráfico e \(x_2\), tempo de cozimento, é o fator de subparcela.

Como as mesmas combinações de tratamento de subparcelas não são executadas em cada parcela inteira, Letsinger et al. (1996) mostram que as estimativas de mínimos quadrados dos parâmetros do modelo, que seriam obtidas pela função rsm, serão insatisfatórias se a razão entre as variâncias do erro do gráfico inteiro e do erro do subparcela for maior que um.

Eles recomendam o uso de estimadores REML, que podem ser calculados na função lmer no pacote lme4. Combinando os coeficientes para os termos lineares, quadráticos e de interação no mesmo vetor, o modelo quadrático geral para um projeto de superfície de resposta completamente aleatório pode ser escrito em forma matricial como \[\begin{equation*} \tag{10.18} \pmb{y}=\pmb{X}\beta+\epsilon\cdot \end{equation*}\]

Quando o planejamento da superfície de resposta é executado como um gráfico dividido (split-plot), como no exemplo da Tabela 10.10, o modelo quadrático geral pode ser escrito na forma \[\begin{equation*} \tag{10.19} \pmb{y}=\pmb{X}\beta+\omega+\epsilon, \end{equation*}\] onde \(\beta\) é o vetor dos coeficientes de regressão para os efeitos da parcela inteira e da subparcela, \(\omega\) é um vetor de erros aleatórios de toda a parcela e \(\epsilon\) é um vetor de erros aleatórios de subparcela.

É assumido que \(\omega+\epsilon\) tem média zero e matriz de variâncias e covariâncias dada por \(\Sigma=\sigma^2\pmb{I}+\sigma_\omega^2\pmb{J}\), onde \(\sigma_\omega^2\) e \(\sigma^2\) são as variâncias das unidades experimentais da parcela inteira e subparcela e \[\begin{equation*} \tag{10.20} \pmb{J}=\begin{pmatrix} \pmb{1}_1\pmb{1}^\top_1 & 0 & \cdots & 0 \\ 0 & \pmb{1}_2\pmb{1}^\top_2 & \cdots & 0 \\ \vdots & \dots & \ddots & \vdots \\ 0 & 0 & \cdots & \pmb{1}_m\pmb{1}^\top_m \end{pmatrix}\cdot \end{equation*}\]

O comprimento de \(\pmb{1}_i\) é \(n_i\), o número de subparcelas executadas dentro da \(i\)-ésima parcela inteira. Os estimadores de mínimos quadrados de \(\beta\) são \(\big(\pmb{X}^\top\pmb{X}\big)^{−1}\pmb{X}^\top\pmb{y}\), enquanto a melhor estimativa linear não viciada de \(\beta\) é \[\begin{equation*} \tag{10.21} \widehat{\beta}=\big(\pmb{X}^\top\Sigma^{-1}\pmb{X}\big)^{-1}\pmb{X}^\top\Sigma^{-1}\pmb{y} \end{equation*}\] e sua matriz de variância e covariância é dada por \[\begin{equation*} \tag{10.22} \mbox{Var}\big(\widehat{\beta}\big)=\big(\pmb{X}^\top\Sigma^{-1}\pmb{X}\big)^{-1}, \end{equation*}\] \(\sigma^2\) e \(\sigma_\omega^2\) devem ser conhecidos para calcular a melhor estimativa linear não viciada de \(\beta\).

No entanto, estimativas de \(\sigma^2\) e \(\sigma_\omega^2\) e das melhores estimativas lineares não viciadas em amostras grandes são encontradas substituindo \(\widehat{\sigma}^2\) e \(\widehat{\sigma}_\omega^2\) na Equação (10.21) podem ser obtidos usando o método REML descrito no Capítulo 5 e disponível no pacote lme4.

Os comandos para ajustar o modelo quadrático geral aos dados da Tabela 10.10 usando a função lmer são mostrados a seguir.

cake.baking = data.frame( cbind(Ovenrun = c(1,1,1,2,2,2,3,3,4,4,4), 
                                x1=c(-1,-1,-1,1,1,1,0,0,0,0,0),
                                x2=c(-1,1,0,-1,1,0,-1,1,0,0,0),
                                y=c(2.7,2.5,2.7,2.9,1.3,2.2,3.7,2.9,2.9,2.8,2.9)))
# REML analysis using lmer
library(lme4)
mmod = lmer(y ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2) + (1|Ovenrun), data = cake.baking)
summary(mmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2) + (1 | Ovenrun)
##    Data: cake.baking
## 
## REML criterion at convergence: -2.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.36458 -0.35677  0.04687  0.47266  0.71354 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ovenrun  (Intercept) 0.1402   0.3745  
##  Residual             0.0025   0.0500  
## Number of obs: 11, groups:  Ovenrun, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  3.13118    0.26666  11.742
## x1          -0.25000    0.26559  -0.941
## x2          -0.43333    0.02041 -21.229
## I(x1^2)     -0.68353    0.37581  -1.819
## I(x2^2)     -0.09648    0.04316  -2.236
## x1:x2       -0.35000    0.02500 -14.000
## 
## Correlation of Fixed Effects:
##         (Intr) x1     x2     I(1^2) I(2^2)
## x1       0.000                            
## x2       0.000  0.000                     
## I(x1^2) -0.703  0.000  0.000              
## I(x2^2) -0.081  0.000  0.000 -0.019       
## x1:x2    0.000  0.000  0.000  0.000  0.000
require(pbkrtest)
# get the KR-approximated degrees of freedom
df.KR = get_Lb_ddf(mmod, fixef(mmod))
# get p-values from the t-distribution using the t-values and approximated
# degrees of freedom
coefs = data.frame(coef(summary(mmod)))
coefs$p.KR = 2 * (1 - pt(abs(coefs$t.value), df.KR))
coefs
##                Estimate Std..Error     t.value       p.KR
## (Intercept)  3.13118490 0.26665801  11.7423244 0.05467867
## x1          -0.25000000 0.26559028  -0.9412995 0.51984868
## x2          -0.43333333 0.02041241 -21.2289111 0.03038050
## I(x1^2)     -0.68352865 0.37580926  -1.8188180 0.32093197
## I(x2^2)     -0.09648438 0.04315832  -2.2355917 0.26872350
## x1:x2       -0.35000000 0.02500000 -14.0000000 0.04593211
# Least Squares Analysis using rsm
library(rsm)
mmodls = rsm(y ~ SO(x1,x2), data=cake.baking)
summary(mmodls)
## 
## Call:
## rsm(formula = y ~ SO(x1, x2), data = cake.baking)
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  2.978947   0.099908 29.8170 7.957e-07 ***
## x1          -0.250000   0.079509 -3.1443  0.025542 *  
## x2          -0.433333   0.079509 -5.4501  0.002826 ** 
## x1:x2       -0.350000   0.097378 -3.5942  0.015638 *  
## x1^2        -0.697368   0.122361 -5.6993  0.002321 ** 
## x2^2         0.152632   0.122361  1.2474  0.267496    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9445, Adjusted R-squared:  0.889 
## F-statistic: 17.01 on 5 and 5 DF,  p-value: 0.003712
## 
## Analysis of Variance Table
## 
## Response: y
##             Df  Sum Sq Mean Sq F value   Pr(>F)
## FO(x1, x2)   2 1.50167 0.75083  19.795 0.004210
## TWI(x1, x2)  1 0.49000 0.49000  12.919 0.015638
## PQ(x1, x2)   2 1.23505 0.61752  16.281 0.006465
## Residuals    5 0.18965 0.03793                 
## Lack of fit  3 0.18298 0.06099  18.298 0.052263
## Pure error   2 0.00667 0.00333                 
## 
## Stationary point of response surface:
##         x1         x2 
## -0.4158277  0.9427722 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1]  0.1872510 -0.7319878
## 
## $vectors
##          [,1]       [,2]
## x1  0.1940643 -0.9809888
## x2 -0.9809888 -0.1940643


Uma comparação das estimativas dos parâmetros obtidas pelo método dos mínimos quadrados e REML é mostrada na Tabela 10.11. Aqui pode ser visto que os efeitos principais lineares estimados são os mesmos para ambos os métodos, mas os erros padrão são muito menores para o efeito de subparcela e maiores para o efeito de parcela inteira usando o método REML mais correto.

As outras estimativas diferem e novamente os testes para os efeitos da subparcela são mais sensíveis usando REML e os testes para os efeitos da parcela inteira são menos sensíveis. As variâncias estimadas do erro REML para toda a parcela e subparcela são mostradas na parte inferior do lado direito da tabela.

Letsinger et al. (1996) mostraram que quando \(\sigma^2_\omega/\sigma^2 < 0.25\) as estimativas de mínimos quadrados estarão razoavelmente próximas das estimativas REML e podem ser usadas na prática. Porém, desde \(\sigma^2_\omega/\sigma^2\) geralmente não é conhecido, o método REML deve ser usado primeiro.

\[ \begin{array}{ccccccc}\hline & \mbox{Mínimos} & \mbox{quadrados} & \mbox{(rsm)} & \mbox{REML} & \mbox{(lmer)} & \\ \mbox{Fator} & \widehat{\beta} & s_\beta & p-\mbox{valor} & \widehat{\beta} & s_\beta & p-\mbox{valor} \\\hline \mbox{(Intercept)} & 2.978947 & 0.099908 & 7.957e-07 & 3.13118490 & 0.26665801 & 0.05467867 \\ x_1 & -0.250000 & 0.079509 & 0.025542 & -0.25000000 & 0.26559028 & 0.51984868 \\ x_2 & -0.433333 & 0.079509 & 0.002826 & -0.43333333 & 0.02041241 & 0.03038050 \\ x_1:x_2 & -0.350000 & 0.097378 & 0.015638 & -0.35000000 & 0.02500000 & 0.04593211 \\ x_1^2 & -0.697368 & 0.122361 & 0.002321 & -0.68352865 & 0.37580926 & 0.32093197 \\ x_2^2 & 0.152632 & 0.122361 & 0.267496 & -0.09648438 & 0.04315832 & 0.26872350 \\\hline & & & & & \widehat{\sigma}_\omega^2= 0.1402 & \widehat{\sigma}^2 = 0.0025 \\\hline \end{array} \] Tabela 10.11: Comparação de mínimos quadrados e estimativas REML para experimento Split-plot de superfície de resposta.

Na Tabela 10.11 mostramos a aproximação de Kenward-Roger para obter os graus de liberdade aproximados e a distribuição \(t\) para obter \(p\)-valores, implementada na função get_Lb_ddf do pacote pbkrtest. Esta função têm como argumentos o modelo ajustado no objeto R mmod e com a função fixef são extraídas as estimativas de efeitos fixos do modelo.

As aproximações de Kenward-Roger produzem \(p\)-valores ligeiramente mais conservadores do que a aproximação normal. No entanto, mesmo neste conjunto de dados não muito grande a aproximação normal é apenas ligeiramente menos conservadora do que esta opção.

Em alguns planejamentos Split-plot, os estimadores de mínimos quadrados são idênticos aos estimadores produzidos pelo REML. Esses delineamentos têm a vantagem de que as estimativas dos parâmetros não dependem dos componentes de variância no delineamento de parcelas subdivididas.

Do ponto de vista prático, seria vantajoso poder utilizar as estimativas de mínimos quadrados em R, uma vez que a função rsm, que utiliza o método de mínimos quadrados, também pode calcular automaticamente a análise canônica, a análise ridge ou produzir valores previstos em uma grade para uso em gráficos de contorno. A função lmer, que produz os estimadores REML, não produz essas análises que são úteis para exploração de superfícies de resposta.

Vining et al. (2005) provaram um teorema de equivalência que mostra que as estimativas de mínimos quadrados de todos os coeficientes de regressão em um projeto Split-plot central de parcela dividida serão as mesmas que as estimativas REML se

  1. o projeto for balanceado, em que cada parcela inteira contém o mesmo número de subparcelas,

  2. os desenhos das subparcelas são ortogonais, embora não necessariamente iguais e

  3. os trechos axiais para as subparcelas são executados em uma única parcela inteira.

Além disso, para obter estimativas separadas das variâncias das unidades experimentais de parcela inteira e subparcela, pelo menos duas parcelas inteiras de pontos centrais em todos os fatores devem ser incluídas. Projetos que possuem essas propriedades serão referidos aqui como “projetos de superfície de resposta Split-plot de estimativas equivalentes” ou EESPRS. Nenhum planejamento composto central da EESPRS é possível quando há apenas um fator de subparcela porque a propriedade (2) não pode ser satisfeita.

Quando há dois ou mais fatores de subparcela, os delineamentos compostos centrais EESPRS podem ser obtidos modificando os delineamentos compostos centrais padrão.

Por exemplo, considere os dois experimentos compostos centrais executados em um arranjo Split-plot na Figura 10.22. Se houver um fator de parcela inteira \(A\), e dois fatores de subparcela \(P\) e \(Q\), o experimento composto central padrão com três pontos centrais mostrados à esquerda pode ser agrupado em cinco parcelas inteiras, classificando no fator de parcela inteira \(A\) conforme exibido na figura.

Figura 10.22: Comparação de Split-plot CCDs.

Como o número de execuções em cada parcela não é igual para o experimento à esquerda, ele não é balanceado e não é um experimento EESPRS. Removendo o ponto central de todo o gráfico que contém os pontos axiais para os fatores de subparcela \(P\) e \(Q\), adicionando réplicas aos gráficos inteiros que contêm os pontos axiais para o fator \(A\) da parcela inteira e adicionando réplicas nas parcelas inteiras que contêm pontos centrais, o plano à direita fica equilibrado.

Como as colunas para os dois fatores de subparcela \(P\) e \(Q\) são ortogonais, ou seja, \(\sum_i p_iq_i = 0\) dentro de cada parcela inteira e todos os pontos axiais para os fatores de subparcela estão na mesma parcela inteira, este é um projeto EESPRS. Parcelas inteiras duplicadas (3 & 5) contendo todos os pontos centrais tornam possível estimar \(\sigma_\omega^2\) e \(\sigma^2\) usando o método REML com a função lmer.

Embora as estimativas de mínimos quadrados e REML dos coeficientes de regressão sejam as mesmas quando um projeto EESPRS é usado, os erros padrão dos coeficientes não serão os mesmos, porque a matriz de covariância das estimativas de mínimos quadrados \(\sigma^2\big(\pmb{X}^\top\pmb{X}\big)^{−1}\), não é igual à covariância das estimativas REML \(\big(\pmb{X}^\top\widehat{\Sigma}\pmb{X}\big)^{−1}\).

Portanto, se for utilizado um projeto EESPRS, a função lmer deve ser usada para fazer testes de hipóteses sobre os parâmetros, a fim de determinar se o modelo pode ser simplificado, mas a função rsm pode ser usada para fazer uma análise canônica, ridge análise ou gráficos de contorno.

Em um projeto EESPRS balanceado, o número de parcelas inteiras é \(m\) e o número de subparcelas dentro de cada parcela inteira é \(n\), perfazendo o número total de execuções \(N = nm\). O número de fatores é \(k = k_1 + k_2\), onde \(k_1\) é o número de fatores do gráfico inteiro e \(k_2\) é o número de fatores do subparcela. Em um planejamento composto central de equivalência de estimativa, o número de execuções em cada parcela inteira será \(2\times k_2\), uma vez que uma parcela inteira deve conter todas as execuções axiais para os fatores da subparcela.

A Tabela 10.12 mostra alguns projetos EESPRS que podem ser construídos modificando projetos compostos centrais padrão. A designação \(k_1(k_2)\) na primeira coluna refere-se ao número de fatores inteiros e subparcelas.

\[ \begin{array}{cccccc}\hline & & & & & \mbox{Fatorial} \\ k_1(k_2) & k & m & n & N & \mbox{Projetos subparcela} \\\hline 1(2) & 3 & 7 & 4 & 28 & 2^2 \\ 1(3) & 4 & 7 & 6 & 42 & 2^{3-1} (\pmb{I}=\pm PQR)+2\mbox{cp} \\ 1(4) & 5 & 7 & 8 & 56 & 2^{4-1} (\pmb{I}=\pm PQRS) \\ 2(2) & 4 & 11 & 4 & 44 & 2^{2} \\ 2(3) & 5 & 11 & 6 & 66 & 2^{3-1} (\pmb{I}=\pm PQR)+2\mbox{cp} \\ 2(4) & 6 & 11 & 8 & 88 & 2^{4-1} (\pmb{I}=\pm PQRS) \\ 3(2) & 5 & 18 & 4 & 72 & 2^{2} \\ 3(3) & 6 & 18 & 6 & 108 & 2^{3-1} (\pmb{I}=\pm PQR)+2\mbox{cp} \\\hline \end{array} \] Tabela 10.12: Projetos de CCD da EESPRS.

O experimento na primeira linha da Tabela 10.11 é o experimento no lado direito da Figura 10.22 que tem um fator de parcela total e dois fatores de subparcela. Neste delineamento, a porção fatorial (parcelas inteiras 2 e 6) possui as quatro combinações de tratamento de um fatorial completo de 22 nos fatores das subparcelas, randomizados para as quatro subparcelas.

Em outros delineamentos mostrados na Tabela 10.11, as subparcelas na parte fatorial do delineamento podem conter um fatorial fracionário nos fatores das subparcelas, ou um fatorial fracionário nos fatores das subparcelas aumentado por dois pontos centrais, conforme indicado no última coluna da tabela.

Usando a Tabela 10.12 como guia, esses projetos podem ser criados nas etapas de programação R. Por exemplo, o código abaixo cria o plano na linha quatro da tabela, que possui dois fatores de gráfico inteiro e dois fatores de subparcela em uma ordem não aleatória.

library("AlgDesign")
# uses gen.factorial function from AlgDesign to create Factorial portion of the design
sp = gen.factorial(2,2,varNames=c("P","Q"))
Wp = gen.factorial(2,2,varNames=c("A","B"))
A = Wp$A
#stacks whole plots
wp = c(rep((1:4),each=4))
A = c(rep((Wp$A),each=4))
B = c(rep((Wp$B),each=4))
Fac = cbind(wp,A,B,rbind(sp,sp,sp,sp))
# Subplot Axial Portion 
A = c(rep(0,4))
B = c(rep(0,4))
P = c(-2,2,0,0)
Q = c(0,0,-2,2)
wp = c(rep(5,4))
spa = cbind(wp,A,B,P,Q)
# Whole Plot Axial Portion
wp = c(rep((6:9),each=4))
A = c(rep(c(-2,2,0,0),each=4))
B = c(rep(c(0,0,-2,2),each=4))
P = c(rep(0,16))
Q = c(rep(0,16))
wpa<-cbind(wp,A,B,P,Q)
# center points
wp = c(rep((10:11),each=4))
A = c(rep(0,8))
B = c(rep(0,8))
C = c(rep(0,8))
D = c(rep(0,8))
wpc = cbind(wp,A,B,P,Q)
SPDs = rbind(Fac,spa,wpa,wpc)
SPDs
##    wp  A  B  P  Q
## 1   1 -1 -1 -1 -1
## 2   1 -1 -1  1 -1
## 3   1 -1 -1 -1  1
## 4   1 -1 -1  1  1
## 5   2  1 -1 -1 -1
## 6   2  1 -1  1 -1
## 7   2  1 -1 -1  1
## 8   2  1 -1  1  1
## 9   3 -1  1 -1 -1
## 10  3 -1  1  1 -1
## 11  3 -1  1 -1  1
## 12  3 -1  1  1  1
## 13  4  1  1 -1 -1
## 14  4  1  1  1 -1
## 15  4  1  1 -1  1
## 16  4  1  1  1  1
## 17  5  0  0 -2  0
## 18  5  0  0  2  0
## 19  5  0  0  0 -2
## 20  5  0  0  0  2
## 21  6 -2  0  0  0
## 22  6 -2  0  0  0
## 23  6 -2  0  0  0
## 24  6 -2  0  0  0
## 25  7  2  0  0  0
## 26  7  2  0  0  0
## 27  7  2  0  0  0
## 28  7  2  0  0  0
## 29  8  0 -2  0  0
## 30  8  0 -2  0  0
## 31  8  0 -2  0  0
## 32  8  0 -2  0  0
## 33  9  0  2  0  0
## 34  9  0  2  0  0
## 35  9  0  2  0  0
## 36  9  0  2  0  0
## 37 10  0  0  0  0
## 38 10  0  0  0  0
## 39 10  0  0  0  0
## 40 10  0  0  0  0
## 41 11  0  0  0  0
## 42 11  0  0  0  0
## 43 11  0  0  0  0
## 44 11  0  0  0  0
## 45 10  0  0  0  0
## 46 10  0  0  0  0
## 47 10  0  0  0  0
## 48 10  0  0  0  0
## 49 11  0  0  0  0
## 50 11  0  0  0  0
## 51 11  0  0  0  0
## 52 11  0  0  0  0


Vining et al. (2005) apresentaram um exemplo de projeto composto central EESPRS com dois fatores de parcela inteira e dois fatores de subparcela. Neste exemplo, um engenheiro estava estudando os efeitos de dois fatores difíceis de alterar \(A\), temperatura do forno da zona 1 e \(B\), temperatura do forno da zona 2 e dois fatores fáceis de alterar \(P\), quantidade de ligante na formulação e \(Q\), velocidade de moagem do lote, na resistência do tubo cerâmico.

Neste exemplo, os pontos axiais estão em \(\pm\) 1, ou seja, planejamento de cubo centrado na face. O dimensionamento em níveis codificados e as medições de resistência resultantes são apresentados na Tabela 10.13. Este projeto usou três parcelas inteiras que consistiam apenas em pontos centrais.

Tabela 10.13: Medidas de projeto e resistência para experimento de tubo cerâmico.

A análise destes dados, deixados para um exercício, mostrará que os coeficientes de regressão para o modelo quadrático geral obtido pelo método dos mínimos quadrados e REML são os mesmos, embora os erros padrão sejam diferentes.

Os projetos Box-Behnken também podem ser modificados para criar projetos EESPRS. Simplesmente classifique o experimento Box-Behnken pelos fatores da parcela inteira e determine o número de subparcelas em cada parcela inteira para combinar o bloco com o número máximo.

Adicione pontos centrais em relação aos fatores de subparcela em parcelas inteiras onde os níveis de fator de parcela inteira não estão no valor central e adicione pelo menos duas parcelas inteiras consistindo inteiramente de pontos centrais. A Tabela 10.14 mostra alguns projetos de EESPRS que podem ser construídos modificando os projetos padrão de Box-Behnken.

\[ \begin{array}{cccccc}\hline k_1(k_2) & k & m & n & N \\\hline 1(2) & 3 & 5 & 4 & 20 \\ 1(3) & 4 & 7 & 6 & 42 \\ 1(4) & 5 & 7 & 8 & 56 \\ 2(2) & 4 & 11 & 4 & 44 \\ 2(3) & 5 & 13 & 6 & 78 \\\hline \end{array} \] Tabela 10.14: Projetos EESPRS BBD.

Usando a Tabela 10.14 como guia, esses projetos também podem ser criados com etapas de programação R. Goos (2006) apontou que muitos dos projetos compostos centrais EESPRS e Box-Behnken são bastante ineficientes devido ao fato de que muitas execuções experimentais são usadas para pontos centrais. Uma vez que a condição para equivalência de estimativa em experimentos Split-plot pode ser escrita em uma equação da forma: \[\begin{equation*} \tag{10.23} \Big(\pmb{I}_n -\pmb{X}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{X}^\top\Big)\pmb{J}\pmb{X} = \pmb{0}_{n\times p}, \end{equation*}\] onde \(\pmb{J}\) é dado na equação 10.20.

Marcharia and Goos (2010) usaram um algoritmo de troca para procurar projetos que fossem equivalentes em estimativa e \(D\)-eficientes. Jones and Goos (2012) melhoraram o algoritmo, implementando uma função de objetivo dupla e encontraram projetos equivalentes de estimativas mais eficientes para muitos dos casos estudados por Marcharia and Goos.

Jones and Goos compilaram um catálogo de 111 planos de superfície de resposta equivalentes a estimativas \(D\)-eficientes produzidos por seu algoritmo. A maioria desses projetos tem três níveis para todos os fatores codificados como −1, 0 e 1, e todos os projetos têm todos os níveis de fator limitados entre −1 e 1.

Este catálogo pode ser recuperado por funções no pacote daewr. A Tabela 10.15 lista as funções para fazer isso.

Tabela 10.15: Funções daewr para recuperar os projetos EESPRS \(D\)-eficientes de Jones and Goos.

Chamar uma das funções sem quaisquer argumentos produz uma tabela listando os nomes dos projetos que podem ser recuperados com essa função, chamar a função com um nome de projeto à medida que o argumento recupera um conjunto de dados contendo o projeto, como mostrado no exemplo a continuação.

library(daewr)
EEw2s3( )
##   
## Catalog of D-efficient Estimation Equivalent RS 
##   Designs for (2 wp factors and  3 sp factors)   
##   
##    Jones and Goos, JQT(2012) pp. 363-374 
##   
## Design Name whole plots sub-plots/whole plot  
## ---------------------------------------- 
## EE24R8WP       8                   3          
## EE28R7WP       7                   4          
## EE32R8WP       8                   4          
## EE35R7WP       7                   5          
## EE40R8WP       8                   5          
## EE42R7WP       7                   6          
## EE48R8WP       8                   6          
##   
## ==> to retrieve a design type EEw2s3('EE24R8WP') etc.
EEw2s3('EE21R7WP')
##    WP w1 w2 s1 s2 s3
## 1   1  1  1 -1 -1  1
## 2   1  1  1  1 -1 -1
## 3   1  1  1 -1  1 -1
## 4   2  0  1  0  1 -1
## 5   2  0  1  1 -1  1
## 6   2  0  1 -1  0  0
## 7   3 -1  0 -1  1  0
## 8   3 -1  0  1 -1 -1
## 9   3 -1  0 -1 -1  1
## 10  4  1 -1  1 -1  1
## 11  4  1 -1 -1  1  1
## 12  4  1 -1  1  1 -1
## 13  5 -1  1 -1 -1 -1
## 14  5 -1  1  1  1  0
## 15  5 -1  1 -1  1  1
## 16  6  1  0  0  0  1
## 17  6  1  0  1  1  1
## 18  6  1  0 -1 -1 -1
## 19  7 -1 -1  0 -1  0
## 20  7 -1 -1 -1  0 -1
## 21  7 -1 -1  1  1  1

10.10 Exercícios