Next: About this document ... Up: Aulas Práticas CE074: Controle Previous: Aula 5: Experimento Fatorial
O objetivo dessa aula é realizar uma análise de superfície de resposta utilizando o R.
Utilizaremos os dados do exemplo estudado em sala de aula para exemplificar uma análise completa.
Inicialmente tem-se os dados referentes ao primeiro experimento, onde o pesquisador iniciou os estudos.
Lembrando que 
representa a variável tempo e 
a variável temperatura. A variável 
representa o rendimento da reação em percentual.
Os dados do experimento podem ser utilizados da seguinte maneira:
x1<-c(-1,-1,1,1,0,0,0,0,0) x2<-c(-1,1,-1,1,0,0,0,0,0) y<-c(39.3,40,40.9,41.5,40.3,40.5,40.7,40.2,40.6) dados1<-data.frame(x1,x2,y)
Em seguida faz-se uma análise de um modelo de regressão para encontrar a relação entre a variável resposta e as variáveis preditoras. Inicialmente, um modelo de primeira ordem é ajustado.
m1.1<-lm(y~x1+x2, data=dados1) summary(m1.1) # estimativa dos efeitos anova(m1.1)
Observe que nos resultados, os componentes inicialmente inseridos no modelo são significativos. Mas, é preciso ainda verificar se podem existir outros componentes como algum efeito de interação ou um algum efeito quadrático.
Call:
lm(formula = y ~ x1 + x2, data = dados1)
Residuals:
Min 1Q Median 3Q Max
-0.244444 -0.044444 0.005556 0.055556 0.255556
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 40.44444 0.05729 705.987 5.45e-16 ***
x1 0.77500 0.08593 9.019 0.000104 ***
x2 0.32500 0.08593 3.782 0.009158 **
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.1719 on 6 degrees of freedom
Multiple R-Squared: 0.941, Adjusted R-squared: 0.9213
F-statistic: 47.82 on 2 and 6 DF, p-value: 0.0002057
O quadro da ANOVA pode ser obtido com o seguinte comando:
> anova(m1.1)
Observe que todos os efeitos são significativos.
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 2.40250 2.40250 81.339 0.0001040 ***
x2 1 0.42250 0.42250 14.304 0.0091581 **
Residuals 6 0.17722 0.02954
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Após a análise de variância do modelo, é necessário verificar se o ajuste do modelo está adequado. Para isso, gera-se um modelo acrescentando-se o componente quadrático e a interação entre os fatores.
m2.1<-lm(y~x1+x2+x1*x2+I(x1^2), data=dados1) anova(m2.1,test="F")
O resultado é
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 2.40250 2.40250 55.8721 0.001713 **
x2 1 0.42250 0.42250 9.8256 0.035030 *
I(x1^2) 1 0.00272 0.00272 0.0633 0.813741 Quadrático
x1:x2 1 0.00250 0.00250 0.0581 0.821316 interação
Residuals 4 0.17200 0.04300
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Tanto o componente quadrático quanto o componente de interação foram não significativos (p>0,05).
Maneira alternativa:
m1.lack<-lm(y~factor(x1)+factor(x2), data=dados1) anova(m1.1,m1.lack,test="F")
O resultado é o seguinte:
Analysis of Variance Table Model 1: y ~ x1 + x2 Model 2: y ~ factor(x1) + factor(x2) Res.Df RSS Df Sum of Sq F Pr(>F) 1 6 0.177222 2 5 0.174500 1 0.002722 0.078 0.7912
Nesse caso, como o p-valor foi não significativo, considera-se que não há falta de ajuste desse modelo.
Ainda, uma outra maneira, é utilizar a função pure.error.anova() do pacote alr3.
Considerando o seguinte modelo:
m2.1<-lm(y~x1+x2+x1*x2, data=dados1)
Pode-se avaliar a falta de ajuste da seguinte maneira:
pure.error.anova(m2.1)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 2.40250 2.40250 55.8721 0.001713 **
x2 1 0.42250 0.42250 9.8256 0.035030 *
x1:x2 1 0.00250 0.00250 0.0581 0.821316
Residuals 5 0.17472 0.03494
Lack of fit 1 0.00272 0.00272 0.0633 0.813741
Pure Error 4 0.17200 0.04300
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Nessa análise deve-se observar o p-valor dos desvios de regressão (Lack of fit). Veja que o p-valor foi não significativo, indicando que não há falta de ajuste nesse modelo.
Tem-se agora, os dados do segundo experimento, onde o pesquisador tem o interesse de ajustar o modelo de segunda ordem:
x1<-c(-1,-1,1,1,0,0,0,0,0) x2<-c(-1,1,-1,1,0,0,0,0,0) y<-c(76.5,77,78,79.5,79.9,80.3,80,79.7,79.8) dados2<-data.frame(x1,x2,y)
A análise do modelo inicial é feita da seguinte maneira
m2.1<-lm(y~x1+x2,data=dados2) summary(m2.1) anova(m2.1)
A análise do ajuste é dada da seguinte forma:
m2.2<-lm(y~x1+x2+x1*x2+I(x1^2), data=dados2) anova(m2.2,test="F")
Observe que nesse caso, o componente quadrático foi significativo (p<0,05).
Usando uma maneira alternativa:
m2.lack<-lm(y~factor(x1)+factor(x2), data=dados2) anova(m2.1,m2.lack,test="F")
Os dados do experimento aumentado são:
x1<-c(-1,-1,1,1,0,0,0,0,0,1.414,-1.414,0,0) x2<-c(-1,1,-1,1,0,0,0,0,0,0,0,1.414,-1.414) y<-c(76.5,77,78,79.5,79.9,80.3,80,79.7,79.8,78.4,75.6,78.5,77) dados3<-data.frame(x1,x2,y)
O modelo pode ser construido da seguinte maneira:
m1.3<-lm(y~x1*x2+I(x1^2)+I(x2^2),data=dados3) summary(m1.3) anova(m1.3)
Observe que os resíduos já são pequenos. O que você acha da falta de ajuste?
Uma das maneiras é utilizando a função persp:
Inicialmente, deve ser gerada uma sequência de valores dos níveis dos fatores.
x<-seq(-1.5,1.5,length=30) y<-seq(-1.5,1.5,length=30)
z<-function(x=x,y=y){y<-m1.3$coef[1]+m1.3$coef[2]*x+m1.3$coef[3]*y+
m1.3$coef[4]*x^2+m1.3$coef[5]*y^2+m1.3$coef[6]*x*y}
#alternativamente, os próprios coeficientes podem ser digitados
#z<-function(x,y){y<-79.94+0.99*x+0.52*y-1.38*x^2-1*y^2+0.25*x*y}
O comando outer prepara os dados para serem utilizados pela função persp
z <- outer(x, y, z) persp(x,y,z,theta = -35, phi = 5, expand = 0.5, col = "gray",xlab="Tempo", ylab="Temperatura",zlab="Rend",scale=T,ticktype = "detailed")
Experimente alterar as opções gráficas!
Pode-se visualizar o comportamento da superfície através de curvas de nível do seguinte modo:
contour(x,y,z,xlab="Tempo",ylab="Temperatura")
Outras representações gráficas podem ser utilizadas.
image(x,y,z,col=topo.colors(12)) image(x,y,z,col=terrain.colors(12)) image(x,y,z,col=heat.colors(12))
Experimente também
image(x,y,z,col=topo.colors(12)) contour(x,y,z,xlab="Tempo",ylab="Temperatura",add=T)
Utilizaremos os procedimentos propostos em sala
b<-matrix(c(m1.3$coef[2],m1.3$coef[3]));b B<-matrix(c(m1.3$coef[4],m1.3$coef[6]/2,m1.3$coef[6]/2,m1.3$coef[5]),ncol=2);B x0<--0.5*solve(B)%*%b ; x0 # coordenadas do ponto máximo # correspondente em tempo e temperatura (var naturais) tempo<-x0[1]*5+85 ;tempo temperatura<-x0[2]*5+175 ; temperatura
Experimente observar o seguinte gráfico:
image(x,y,z,col=topo.colors(12)) contour(x,y,z,xlab="Tempo",ylab="Temperatura",add=T) points(x0[1],x0[2],col="red")
O ponto "máximo" é dado por
y0<-m1.3$coef[1]+1/2*t(x0)%*%b
Para analisar o ponto estacionário devemos observar os sinais das raízes.
raizes<-eigen(B)
> raizes
$values
[1] -0.9634986 -1.4142867
$vectors
[,1] [,2]
[1,] 0.2897174 0.9571122
[2,] 0.9571122 -0.2897174
Como os sinais são negativos, o ponto é de máximo.