Adição de pontos centrais a um planejamento \(2^k\)
Nos experimentos \(2^k\), pelo fato de usarmos apenas 2 níveis de cada fator, estamos automaticamente assumindo que a relação de dependência da resposta para os fatores é linear. Existe então uma preocupação sobre essa suposição de linearidade.
Podemos notar que, ao se adicionar um termo de interação à um modelo de efeitos principais de primeira ordem, resultando em
- \[ y = \beta_0 + \sum_{j=1}^{k} \beta_{j} x_j + \mathop{\sum\sum}\limits_{i < j} \beta_{ij} x_i x_j + \epsilon \]
então temos um modelo capaz de representar uma curvatura na variável resposta. Essa curvatura é resultado da “torção” no plano induzida pelos termos da interação \(\beta_{ij} x_i x_j\).
No entanto, podem existir situações onde a curvatura na variável resposta não é modelada adequadamente pela equação 1 acima. Nestes casos, devemos considerar um modelo com
- \[ y = \beta_0 + \sum_{j=1}^{k} \beta_{j} x_j + \mathop{\sum\sum}\limits_{i < j} \beta_{ij} x_i x_j + \sum_{j=1}^{k} \beta_{jj} x_{j}^{2} + \epsilon \]
onde \(\beta_{jj}\) representa efeitos quadráticos ou de segunda ordem. Este modelo também é chamado de superfície de resposta de segunda ordem.
Uma forma de testar se existe curvatura (“proteção contra curvatura”) é através da adição de pontos centrais ao planejamento \(2^k\). Como o próprio nome sugere, esses pontos serão os valores (0, 0) na escala codificada dos níveis (-1, 1). Portanto, o processo consiste em realizar \(n_c\) réplicas no ponto \(x_i\) (\(i = 1, 2, \ldots, k\)). Uma razão importante para adicionar as corridas replicadas no centro do planejamento, é que pontos centrais não interferem nas estimativas dos efeitos de planejamentos \(2^k\). Quando adicionamos os pontos centrais, consideramos que os \(k\) fatores são quantitativos.
Para ilustrar a abordagem, considere um experimento \(2^2\) com uma observação em cada um dos pontos fatoriais: (-,-), (+,-), (-,+), (+,+), e \(n_c\) observações nos pontos centrais (0,0). Seja \(\bar{y}_F\) a média das quatro corridas nos quatro pontos fatoriais, e \(\bar{y}_C\) a média das \(n_C\) corridas no ponto central. Se a diferença \(\bar{y}_F - \bar{y}_C\) for pequena, os pontos centrais estarão no plano ou próximo do plano passando através dos pontos fatoriais, não havendo portanto curvatura. Por outro lado, se \(\bar{y}_F - \bar{y}_C\) for grande, então há indícios de que existe curvatura. A soma de quadrados da curvatura, com um grau de liberdade é dada por
- \[ \begin{align} SQ_{curv} &= \left( \frac{\bar{y}_F - \bar{y}_C}{\sqrt{\frac{1}{n_F} + \frac{1}{n_C}}} \right)^2 \\ \\ &= \frac{n_F n_C (\bar{y}_F - \bar{y}_C)^2}{n_F + n_C} \end{align} \]
onde \(n_F\) é o número de pontos do planejamento fatorial. Essa soma de quadrados deve ser comparada com a média quadrática do erro (como existe um grau de liberdade, então \(SQ_{curv} = MQ_{curv}\)) para testar a curvatura. Note que quando a equação 3) acima for dividida por \(\hat{\sigma}^2 = MQ_{Res}\), o resultado será similar ao quadrado da estatística \(t\) usada para comparar duas médias.
Mais especificamente, quando pontos centrais são adicionados à um experimento \(2^k\), então o teste para curvatura (usando a equação acima) realmente testa as hipóteses
\[ \begin{align} \text{H}_0 &: \sum_{j=1}^{k} \beta_{jj} = 0 \\ \text{H}_1 &: \sum_{j=1}^{k} \beta_{jj} \neq 0 \end{align} \]
Além disso, se os pontos fatoriais do planejamento não forem replicados, poderemos usar os \(n_C\) pontos centrais para construir uma estimativa de erro com \(n_C - 1\) graus de liberdade.
Com a adição de pontos centrais é possível ter uma estimativa de erro (variância) pura, ou seja, que independe do modelo especificado. A idéia básica é decompor a soma de quadrados dos resíduos em soma de quadrado pura e soma de quadrados devido à falta de ajuste.
O teste da falta de ajuste é importante para situar-se na região experimental. No caso de não haver falta de ajuste entende-se que o modelo de primeira ordem (1) que estima efeitos principais e interações apenas, é adequado para descrever o fenômeno. Por outro lado, a existência de falta de ajuste indica que para pelo menos um dos \(k\) fatores existe uma relação não linear e que este fator está perto da sua região experimental ótima.
Quando isso acontecer, os próximos experimentos devem ser planejados de forma a pertimitir estimação de termos quadráticos para capturar essa não linearidade da relação \(y \sim x_1 + x_2 + \cdots + x_k\) bem como permitir encontrar pontos de ótimo da relação, que é o nosso objetivo principal.
Sobre a falta de ajuste em modelos lineares
De maneira geral, para verificar se o modelo linear é adequado para descrever a relação entre \(X\) e \(Y\), devemos realizar um teste para falta de ajuste.
A falta de ajuste pode ser verificada quando existem medições repetidas de \(Y\) para um certo nível de \(X\). Por exemplo, se \(X_i\) for considerada em uma escala discreta com \(k\) níveis (\(i = 1, 2, \ldots, k\)), podemos observar \(n_i\) repetições de \(Y_{ij}\) (\(j = 1, 2, \ldots, n_i\)).
No caso específico de um experimento fatorial \(2^2\), por exemplo, temos o fator \(A\) com 2 níveis, \(B\) com 2 níveis, onde cada nível (\(-1\) e \(1\)) são os pontos na escala discreta. Considere que não há repetições, então existem apenas quatro observações deste experimento fatorial. Com a adição de um ponto central (\(0,0\)), geralmente temos \(n_C\) repetições, que seriam as medidas repetidas para este ponto.
Quando existem repetições, é possível calcular o erro puro, ou seja, o desvio de cada observação \(Y_{ij}\), com relação à média do nível \(i\), \(\bar{y}_i\).
Novamente fazendo associação com um experimento fatorial com ponto central, nós podemos usar as medidads repetidas do ponto central para estimar o erro experimental. Note que esta é uma grande vantagem da adição de pontos centrais, pois mesmo em experimentos fatorias onde não há repetições para os pontos dos fatores, se houver medidas repetidas apenas no ponto central, podemos ter uma estimativa do erro.
Note que o erro puro não depende de nenhum modelo, pois mede apenas os desvios das observações em relação à média. Por esse motivo, ele pode ser comparado com os desvios de um modelo de regressão para verificação da qualidade do ajuste. A ideia é que se o modelo linear é adequado, então os desvios das observações em relação à média (erro puro) devem ser próximos dos desvios das observações em relação ao modelo ajustado.
Se não houverem repetições nos níveis de \(X_i\) (ou no ponto central), então a verificação da suposição de linearidade deve ser verificada através dos resíduos do modleo apenas. (Mesmo com a realização de um teste de falta de ajuste, os resíduos sempre devem ser verificados).
Para representar o erro puro e a falta de ajuste, podemos realizar a decomposição da soma de quadrados dos resíduos, quando existem repetições para \(X_i\), em
\[ \begin{align} A &= B + C \\ (y_{ij} - \hat{y}_i) &= (y_{ij} - \bar{y}_i) + (\bar{y}_i - \hat{y}_i) \\ \text{Resíduo} &= \text{Erro puro} + \text{Desvio de regressão (falta de ajuste)} \end{align} \]
Fazendo a soma de quadrados destes termos, temos
\[ \begin{align} \sum_{i=1}^{k}\sum_{j=1}^{n_i} (y_{ij} - \hat{y}_i)^2 &= \sum_{i=1}^{k}\sum_{j=1}^{n_i} (y_{ij} - \bar{y}_i)^2 + \sum_{i=1}^{k}\sum_{j=1}^{n_i} (\bar{y}_i - \hat{y}_i)^2 \\ &= \sum_{i=1}^{k}\sum_{j=1}^{n_i} (y_{ij} - \bar{y}_i)^2 + \sum_{i=1}^{k} n_i (\bar{y}_i - \hat{y}_i)^2 \end{align} \]
Com isso, podemos ver que um modelo “bom” (onde a suposição de linearidade é adequada), será aquele em que a falta de ajuste (ou desvio de regeressão), será próxima de zero, pois em um modleo correto, esperamos que
\[ \hat{y}_i = \bar{y}_i \]
Com isso, se o modelo for correto, então \(QMRes\) provê uma estimativa não viesada para \(\sigma^2\). Se o modelo for incorreto, então \(QMRes\) será uma estimativa viesada para \(\sigma^2\), pois haverá um desvio sistemático (viés) entre o modelo correto (desconhecido), e o modelo incorreto utilizado:
\[ \sigma^2 + \text{B}_i \]
onde
\[ \text{B}_i = \text{E}(Y_i) - \text{E}(\hat{Y}_i) \]
é o viés. Usar o modelo correto, levará, portanto, à estimativas viciadas das variâncias dos estimadores, e por consequência, os intervalos de confiança e os testes de hipótese serão incorretos.
Faz sentido então, testarmos as hipóteses
\[ \begin{align} \text{H}_0 &: \text{o modelo é adequado (não há falta de ajuste)} \\ \text{H}_1 &: \text{o modelo não é adequado (há falta de ajuste)} \end{align} \]
com base na \(SQep\) (soma de quadrados do erro puro) e \(SQlof\) (soma de quadrados da falta de ajuste). Como já vimos, a \(SQep\) é calculada como
\[ SQep = \sum_{i=1}^{k}\sum_{j=1}^{n_i} (y_{ij} - \bar{y}_i)^2 \]
com \((n-k)\) graus de liberdade, pois são necessárias \(k\) médias para o cálculo. No caso do experimento fatorial \(2^2\) com \(n_C\) pontos centrais, esta soma de quadrados é
\[ SQep = \sum_{i=1}^{n_C} (y_{i} - \bar{y}_C)^2 \]
onde \(\bar{y}_C\) é a média dos pontos centrais. Como aqui só há uma média, então existem \(n_C - 1\) graus de liberdade.
Como vimos na equação (3) acima, a soma de quadrados da curvatura, ou da falta de ajuste, é
\[ SQ_{curv} = \frac{n_F n_C (\bar{y}_F - \bar{y}_C)^2}{n_F + n_C} \]
com um grau de liberdade. Com isso, a hipótese pode ser testada com
\[ \frac{SQcurv/1}{SQep/(n_C - 1)} = \frac{MQcurv}{MQep} \]
que, sob \(H_0\), possui distribuição \(F\) com \(1\), e \(n_C - 1\) graus de liberdade.
Exemplos
##----------------------------------------------------------------------
## Análise do 2²+5 pontos centrais
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/montgomery_ex_6-6.txt"
dados <- read.csv(url, sep = "")
## Codifica A e B
dados <- transform(dados,
Ac = (A - mean(A))/5,
Bc = (B - mean(B))/5)
## Definições do experimento
k <- 2
n <- 1
##----------------------------------------------------------------------
## Montando a tabela de sinais
tab <- model.matrix(~ Ac * Bc, data = dados)
tab (Intercept) Ac Bc Ac:Bc
1 1 -1 -1 1
2 1 1 -1 -1
3 1 -1 1 -1
4 1 1 1 1
5 1 0 0 0
6 1 0 0 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
attr(,"assign")
[1] 0 1 2 3##----------------------------------------------------------------------
## Estimativa dos efeitos
## Podemos usar as colunas da tabela de sinais para calcular os
## contrastes
y <- dados$y
## Contrastes
(contr <- t(tab[, -1]) %*% y) [,1]
Ac 3.1
Bc 1.3
Ac:Bc -0.1## Efeitos = contraste/(n2^{k-1})
(ef <- contr/(n * 2^(k - 1))) [,1]
Ac 1.55
Bc 0.65
Ac:Bc -0.05##----------------------------------------------------------------------
## Somas de quadrados, Quadrados médios e teste F
## SQ fatores
(SQf <- contr^2/(n * 2^k)) [,1]
Ac 2.4025
Bc 0.4225
Ac:Bc 0.0025## Separa os pontos fatoriais dos centrais
yf <- dados$y[dados$Ac != 0] # y dos pontos fatoriais
yc <- dados$y[dados$Ac == 0] # y do ponto central (0,0)
nf <- length(yf)
nc <- length(yc)
## SQ da curvatura (falta de ajuste)
## lof = lack of fit
(SQlof <- (nf * nc * (mean(yc) - mean(yf))^2)/(nf + nc))[1] 0.002722222## SQ do erro puro
(SQpuro <- sum((yc - mean(yc))^2))[1] 0.172## Teste F
(Flof <- SQlof/(SQpuro/(nc - 1)))[1] 0.06330749pf(Flof, 1, nc - 1, lower.tail = FALSE)[1] 0.8137408##======================================================================
## Como fazer pela lm()
## Contraste que define a falta de ajuste.
dados$lof <- dados$Ac^2
## Ajuste do modelo
## Modleo simples, assume linearidade
m0 <- lm(y ~ Ac + Bc + Ac:Bc, data = dados)
anova(m0)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
Ac 1 2.40250 2.40250 68.7520 0.0004166 ***
Bc 1 0.42250 0.42250 12.0906 0.0177127 *
Ac:Bc 1 0.00250 0.00250 0.0715 0.7997870
Residuals 5 0.17472 0.03494
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Modelo ocm um contraste para lof
m1 <- lm(y ~ Ac + Bc + lof + Ac:Bc, data = dados)
anova(m1)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
Ac 1 2.40250 2.40250 55.8721 0.001713 **
Bc 1 0.42250 0.42250 9.8256 0.035030 *
lof 1 0.00272 0.00272 0.0633 0.813741
Ac:Bc 1 0.00250 0.00250 0.0581 0.821316
Residuals 4 0.17200 0.04300
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## O ideal é que a SQlof seja a última linha do quadro de anova, porém,
## como nesses modelos existe ortogonalidade entre os efeitos, isso deixa
## de ter importância.
## TRV para testar se lof é significativa
anova(m0, m1)Analysis of Variance Table
Model 1: y ~ Ac + Bc + Ac:Bc
Model 2: y ~ Ac + Bc + lof + Ac:Bc
Res.Df RSS Df Sum of Sq F Pr(>F)
1 5 0.17472
2 4 0.17200 1 0.0027222 0.0633 0.8137##======================================================================
## Exercício 14-23, pg 356, Montgomery, Estatística aplicada e
## probabilidade para Engenheiros, 4 ed. Considere um experimento
## fatorial 2² com 4 pontos centrais. Os dados são (1)=21, a=125, b=154,
## ab=352 e as respostas no ponto central são 92, 130, 98, 152. Calcule
## uma ANOVA com a soma de quadrados para a curvatura (lof) e faça o
## teste F, use alpha=0.05.
da <- data.frame(A = c(-1,1,-1,1,0,0,0,0),
B = c(-1,-1,1,1,0,0,0,0),
y = c(21,125,154,352,92,130,98,152))
## Tabela de sinais
X <- model.matrix(~A*B, data=da)
X (Intercept) A B A:B
1 1 -1 -1 1
2 1 1 -1 -1
3 1 -1 1 -1
4 1 1 1 1
5 1 0 0 0
6 1 0 0 0
7 1 0 0 0
8 1 0 0 0
attr(,"assign")
[1] 0 1 2 3nr <- 1; k <- 2 # nº de repetições e fatores
contr <- t(X[,-1])%*%da$y # contrastes
SQf <- contr^2/(n * 2^k) # soma de quadrado dos fatores
SQf [,1]
A 22801
B 32400
A:B 2209yf <- da$y[da$A!=0] # y dos pontos fatoriais
yc <- da$y[da$A==0] # y do ponto central (0,0)
nf <- length(yf)
nc <- length(yc)
SQlof <- nf*nc*(mean(yc)-mean(yf))^2/(nf+nc) # soma de quadrado devido à lof
SQpura <- sum((yc-mean(yc))^2) # soma de quadrados do erro puro
Flof <- SQlof/(SQpura/(nc-1))
pf(Flof, 1, nc-1, lower=FALSE) # não significativo[1] 0.1087917Ff <- SQf/(SQpura/(nc-1))
pf(Ff, 1, nc-1, lower=FALSE) [,1]
A 0.012671294
B 0.007740781
A:B 0.193497730anova(lm(y ~ A*B, data = da))Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
A 1 22801 22801 14.193 0.01965 *
B 1 32400 32400 20.168 0.01090 *
A:B 1 2209 2209 1.375 0.30602
Residuals 4 6426 1606
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1