La superficie de respuesta (RSM) es un tipo de diseño experimental que permite obtener información para interpolar los resultados y calcular una superficie que permite aproximar los valores del resultado en un rango de los valores de respuesta. RSM se aplica cuando hay evidencia o duda razonable de que haya relaciones de segundo orden significativas.
Diseño experimental
El DoE se hace con base en la librería {rsm}.
library (rsm)
trsm <- ccd(basis = 2, randomize = F,n0 = c(3,3),
coding = list(x1~(RA-163)/15, x2~(CE-250)/35))
Resultados del experimento
Los resultados se agregan como una columna adicional al data.frame y se recodifica la tabla.
Y7 <- c(199, 370, 292, 549, 356, 352, 357, 195, 486, 258, 450, 364, 369, 366)
trsm <- data.frame(trsm,Y7)
trsm <- coded.data(trsm, x1~(ra-163)/15,
x2~(ce-250)/35)
trsm
run.order std.order ra ce Block Y7
1 1 1 148.0000 215.0000 1 199
2 2 2 178.0000 215.0000 1 370
3 3 3 148.0000 285.0000 1 292
4 4 4 178.0000 285.0000 1 549
5 5 5 163.0000 250.0000 1 356
6 6 6 163.0000 250.0000 1 352
7 7 7 163.0000 250.0000 1 357
8 1 1 141.7868 250.0000 2 195
9 2 2 184.2132 250.0000 2 486
10 3 3 163.0000 200.5025 2 258
11 4 4 163.0000 299.4975 2 450
12 5 5 163.0000 250.0000 2 364
13 6 6 163.0000 250.0000 2 369
14 7 7 163.0000 250.0000 2 366
Data are stored in coded form using these coding formulas ...
x1 ~ (ra - 163)/15
x2 ~ (ce - 250)/35
Análisis de los resultados.
first <- rsm(Y7~SO(x1,x2),data = trsm)
summary(first)
Call:
rsm(formula = Y7 ~ SO(x1, x2), data = trsm)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 360.6667 2.5346 142.2957 6.654e-15 ***
x1 104.9420 2.1951 47.8084 4.052e-11 ***
x2 67.9411 2.1951 30.9519 1.290e-09 ***
x1:x2 21.5000 3.1043 6.9259 0.0001213 ***
x1^2 -8.7708 2.2847 -3.8390 0.0049544 **
x2^2 -2.0208 2.2847 -0.8845 0.4022195
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.9976, Adjusted R-squared: 0.9961
F-statistic: 661.3 on 5 and 8 DF, p-value: 3.053e-10
Analysis of Variance Table
Response: Y7
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2) 2 125031 62515 1621.8342 3.664e-11
TWI(x1, x2) 1 1849 1849 47.9686 0.0001213
PQ(x1, x2) 2 582 291 7.5435 0.0144177
Residuals 8 308 39
Lack of fit 3 89 30 0.6766 0.6027711
Pure error 5 219 44
Stationary point of response surface:
x1 x2
-4.816311 -8.810611
Stationary point in original units:
ra ce
90.75534 -58.37138
Eigenanalysis:
eigen() decomposition
$values
[1] 5.871514 -16.663181
$vectors
[,1] [,2]
x1 -0.5918031 -0.8060825
x2 -0.8060825 0.5918031
best <- rsm(Y7~FO(x1,x2)+TWI(x1,x2)+PQ(x1),trsm)
summary(best)
Call:
rsm(formula = Y7 ~ FO(x1, x2) + TWI(x1, x2) + PQ(x1), data = trsm)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 359.4231 2.0833 172.5272 < 2.2e-16 ***
x1 104.9420 2.1684 48.3972 3.441e-12 ***
x2 67.9411 2.1684 31.3331 1.685e-10 ***
x1:x2 21.5000 3.0665 7.0112 6.247e-05 ***
x1^2 -8.6154 2.2502 -3.8287 0.004036 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.9974, Adjusted R-squared: 0.9962
F-statistic: 847 on 4 and 9 DF, p-value: 1.392e-11
Analysis of Variance Table
Response: Y7
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2) 2 125031 62515 1662.0248 2.763e-12
TWI(x1, x2) 1 1849 1849 49.1573 6.247e-05
PQ(x1) 1 551 551 14.6591 0.004036
Residuals 9 339 38
Lack of fit 4 119 30 0.6793 0.635323
Pure error 5 219 44
Stationary point of response surface:
x1 x2
-3.160052 -7.413588
Stationary point in original units:
ra ce
115.599215 -9.475595
Eigenanalysis:
eigen() decomposition
$values
[1] 7.273271 -15.888656
$vectors
[,1] [,2]
x1 -0.5603734 -0.8282401
x2 -0.8282401 0.5603734
contour (best, ~x1+x2, image = T)
persp (best, ~x1+x2, theta = -15)
persp (best, ~x1+x2, theta = 15, contours = "col", col = rainbow(40))
persp (best, ~x1+x2, theta = 30, contours = "col", col = rainbow(40))
NA
NA
Podemos usar la gráfica de contornos para determinar valores de los predictores que produzcan valores deseados de la respuesta.
contour (best, ~x1+x2, image = T)
# 250
point <- data.frame(ra=c(140:144),
ce=c(rep(300,5)))
Yp <- predict(best, coded.data(point,x1~(ra-163)/15,
x2~(ce-250)/35))
Param250 <- data.frame(point,Yp)
Param250
# 400
point <- data.frame(ra=c(156:160),
ce=c(rep(300,5)))
Yp <- predict(best, coded.data(point,x1~(ra-163)/15,
x2~(ce-250)/35))
Param400 <- data.frame(point,Yp)
Param400
# 550
point <- data.frame(ra=c(173:177),
ce=c(rep(300,5)))
Yp <- predict(best, coded.data(point,x1~(ra-163)/15,
x2~(ce-250)/35))
Param550 <- data.frame(point,Yp)
Param550
Ejercicio 1: Determinar un valor de ce para obtener una distancia de 550, con un ra=180
point <- data.frame(ra=c(rep(180,5)),
ce=c(280:284))
Yp <- predict(best, coded.data(point,x1~(ra-163)/15,
x2~(ce-250)/35))
Param550 <- data.frame(point,Yp)
Param550
Ejercicio 2: Determina un valor de cup elevation (ce) que con un valor de tu elección de ra produzcan un disparo con distancia de 500.
library (rsm)
heli <- ccd(basis = 2, randomize = F,n0 = c(3,3),
coding = list(x1~(angle-10)/2, x2~(wings-75)/2))
heli
run.order std.order angle wings Block
1 1 1 8.000000 73.00000 1
2 2 2 12.000000 73.00000 1
3 3 3 8.000000 77.00000 1
4 4 4 12.000000 77.00000 1
5 5 5 10.000000 75.00000 1
6 6 6 10.000000 75.00000 1
7 7 7 10.000000 75.00000 1
8 1 1 7.171573 75.00000 2
9 2 2 12.828427 75.00000 2
10 3 3 10.000000 72.17157 2
11 4 4 10.000000 77.82843 2
12 5 5 10.000000 75.00000 2
13 6 6 10.000000 75.00000 2
14 7 7 10.000000 75.00000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (angle - 10)/2
x2 ~ (wings - 75)/2
Se llevará a cabo un experimento para evaluar el impacto de la temperatura y la concentración en el yield de en un reactor químico. La temperatura se explorará de 50 a 100 grados centrígrados. La concentración se probará de 25% a 50%. Se sospecha una respuesta de segundo orden. Considerando que es posible operar el reactor fuera de los parámetros indicados, recomiende los valores no codificados para las corridas experimentales necesarias.
library (rsm)
reactor <- ccd(basis = 2, randomize = F, n0 = c(3,3),
coding = list(x1~(temp-75)/25,x2~(conc-37.5)/12.5))
reactor
run.order std.order temp conc Block
1 1 1 50.00000 25.00000 1
2 2 2 100.00000 25.00000 1
3 3 3 50.00000 50.00000 1
4 4 4 100.00000 50.00000 1
5 5 5 75.00000 37.50000 1
6 6 6 75.00000 37.50000 1
7 7 7 75.00000 37.50000 1
8 1 1 39.64466 37.50000 2
9 2 2 110.35534 37.50000 2
10 3 3 75.00000 19.82233 2
11 4 4 75.00000 55.17767 2
12 5 5 75.00000 37.50000 2
13 6 6 75.00000 37.50000 2
14 7 7 75.00000 37.50000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (conc - 37.5)/12.5
library (rsm)
reactor.i <- ccd(basis = 2, randomize = F, n0 = c(3,3), inscribed = T,
coding = list(x1~(temp-75)/25,x2~(conc-37.5)/12.5))
reactor.i
run.order std.order temp conc Block
1 1 1 57.32233 28.66117 1
2 2 2 92.67767 28.66117 1
3 3 3 57.32233 46.33883 1
4 4 4 92.67767 46.33883 1
5 5 5 75.00000 37.50000 1
6 6 6 75.00000 37.50000 1
7 7 7 75.00000 37.50000 1
8 1 1 50.00000 37.50000 2
9 2 2 100.00000 37.50000 2
10 3 3 75.00000 25.00000 2
11 4 4 75.00000 50.00000 2
12 5 5 75.00000 37.50000 2
13 6 6 75.00000 37.50000 2
14 7 7 75.00000 37.50000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (conc - 37.5)/12.5
library (rsm)
reactor.f <- ccd(basis = 2, randomize = F, n0 = c(3,3), alpha="faces",
coding = list(x1~(temp-75)/25,x2~(conc-37.5)/12.5))
reactor.f
run.order std.order temp conc Block
1 1 1 50 25.0 1
2 2 2 100 25.0 1
3 3 3 50 50.0 1
4 4 4 100 50.0 1
5 5 5 75 37.5 1
6 6 6 75 37.5 1
7 7 7 75 37.5 1
8 1 1 50 37.5 2
9 2 2 100 37.5 2
10 3 3 75 25.0 2
11 4 4 75 50.0 2
12 5 5 75 37.5 2
13 6 6 75 37.5 2
14 7 7 75 37.5 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (conc - 37.5)/12.5
Con el helicóptero JPL implementa un CCC, CCI y CCF
library (rsm)
heli <- ccd(basis = 2, randomize = F, n0 = c(3,3),
coding = list(x1~(la-105)/2,x2~(aa-10)/2))
heli
run.order std.order la aa Block
1 1 1 103.0000 8.000000 1
2 2 2 107.0000 8.000000 1
3 3 3 103.0000 12.000000 1
4 4 4 107.0000 12.000000 1
5 5 5 105.0000 10.000000 1
6 6 6 105.0000 10.000000 1
7 7 7 105.0000 10.000000 1
8 1 1 102.1716 10.000000 2
9 2 2 107.8284 10.000000 2
10 3 3 105.0000 7.171573 2
11 4 4 105.0000 12.828427 2
12 5 5 105.0000 10.000000 2
13 6 6 105.0000 10.000000 2
14 7 7 105.0000 10.000000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
library (rsm)
heli.i <- ccd(basis = 2, randomize = F, n0 = c(3,3), inscribed = T,
coding = list(x1~(la-105)/2,x2~(aa-10)/2))
heli.i
run.order std.order la aa Block
1 1 1 103.5858 8.585786 1
2 2 2 106.4142 8.585786 1
3 3 3 103.5858 11.414214 1
4 4 4 106.4142 11.414214 1
5 5 5 105.0000 10.000000 1
6 6 6 105.0000 10.000000 1
7 7 7 105.0000 10.000000 1
8 1 1 103.0000 10.000000 2
9 2 2 107.0000 10.000000 2
10 3 3 105.0000 8.000000 2
11 4 4 105.0000 12.000000 2
12 5 5 105.0000 10.000000 2
13 6 6 105.0000 10.000000 2
14 7 7 105.0000 10.000000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
library (rsm)
heli.f <- ccd(basis = 2, randomize = F, n0 = c(3,3), alpha="faces",
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
heli.f
run.order std.order la aa Block
1 1 1 103 8 1
2 2 2 107 8 1
3 3 3 103 12 1
4 4 4 107 12 1
5 5 5 105 10 1
6 6 6 105 10 1
7 7 7 105 10 1
8 1 1 103 10 2
9 2 2 107 10 2
10 3 3 105 8 2
11 4 4 105 12 2
12 5 5 105 10 2
13 6 6 105 10 2
14 7 7 105 10 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
##CCC
library (rsm)
rq.ccc <- ccd(basis = 2, randomize = F, n0 = c(3,3),
coding = list(x1~(temp-75)/25,x2~(con-37.5)/12.5))
rq.ccc
run.order std.order temp con Block
1 1 1 50.00000 25.00000 1
2 2 2 100.00000 25.00000 1
3 3 3 50.00000 50.00000 1
4 4 4 100.00000 50.00000 1
5 5 5 75.00000 37.50000 1
6 6 6 75.00000 37.50000 1
7 7 7 75.00000 37.50000 1
8 1 1 39.64466 37.50000 2
9 2 2 110.35534 37.50000 2
10 3 3 75.00000 19.82233 2
11 4 4 75.00000 55.17767 2
12 5 5 75.00000 37.50000 2
13 6 6 75.00000 37.50000 2
14 7 7 75.00000 37.50000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (con - 37.5)/12.5
##CCI
library (rsm)
rq.cci <- ccd(basis = 2, randomize = F, n0 = c(3,3), inscribed = T,
coding = list(x1~(temp-75)/25,x2~(con-37.5)/12.5))
rq.cci
run.order std.order temp con Block
1 1 1 57.32233 28.66117 1
2 2 2 92.67767 28.66117 1
3 3 3 57.32233 46.33883 1
4 4 4 92.67767 46.33883 1
5 5 5 75.00000 37.50000 1
6 6 6 75.00000 37.50000 1
7 7 7 75.00000 37.50000 1
8 1 1 50.00000 37.50000 2
9 2 2 100.00000 37.50000 2
10 3 3 75.00000 25.00000 2
11 4 4 75.00000 50.00000 2
12 5 5 75.00000 37.50000 2
13 6 6 75.00000 37.50000 2
14 7 7 75.00000 37.50000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (con - 37.5)/12.5
##CCF
library (rsm)
rq.ccf <- ccd(basis = 2, randomize = F, n0 = c(3,3), alpha="faces",
coding = list(x1~(temp-75)/25,x2~(con-37.5)/12.5))
rq.ccf
run.order std.order temp con Block
1 1 1 50 25.0 1
2 2 2 100 25.0 1
3 3 3 50 50.0 1
4 4 4 100 50.0 1
5 5 5 75 37.5 1
6 6 6 75 37.5 1
7 7 7 75 37.5 1
8 1 1 50 37.5 2
9 2 2 100 37.5 2
10 3 3 75 25.0 2
11 4 4 75 50.0 2
12 5 5 75 37.5 2
13 6 6 75 37.5 2
14 7 7 75 37.5 2
Data are stored in coded form using these coding formulas ...
x1 ~ (temp - 75)/25
x2 ~ (con - 37.5)/12.5
Demostración de aleatorización de tabla experimental
library (rsm)
heli.f <- ccd(basis = 2, n0 = c(1,1), alpha="faces",
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
heli.f
run.order std.order la aa Block
1 1 4 107 12 1
2 2 5 105 10 1
3 3 1 103 8 1
4 4 3 103 12 1
5 5 2 107 8 1
6 1 1 103 10 2
7 2 3 105 8 2
8 3 5 105 10 2
9 4 2 107 10 2
10 5 4 105 12 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
2da corrida
library (rsm)
heli.1 <- ccd(basis = 2, n0 = c(1,1), alpha="faces",
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
heli.1
run.order std.order la aa Block
1 1 5 105 10 1
2 2 3 103 12 1
3 3 4 107 12 1
4 4 2 107 8 1
5 5 1 103 8 1
6 1 1 103 10 2
7 2 5 105 10 2
8 3 2 107 10 2
9 4 4 105 12 2
10 5 3 105 8 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
#CCF
library (rsm)
h.ccf <- ccd(basis = 2, randomize = F, n0 = c(2,2), alpha="faces",
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
h.ccf
run.order std.order la aa Block
1 1 1 103 8 1
2 2 2 107 8 1
3 3 3 103 12 1
4 4 4 107 12 1
5 5 5 105 10 1
6 6 6 105 10 1
7 1 1 103 10 2
8 2 2 107 10 2
9 3 3 105 8 2
10 4 4 105 12 2
11 5 5 105 10 2
12 6 6 105 10 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
#CCC
h.ccc <- ccd(basis = 2, randomize = F, n0 = c(2,2),
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
h.ccc
run.order std.order la aa Block
1 1 1 103.0000 8.000000 1
2 2 2 107.0000 8.000000 1
3 3 3 103.0000 12.000000 1
4 4 4 107.0000 12.000000 1
5 5 5 105.0000 10.000000 1
6 6 6 105.0000 10.000000 1
7 1 1 102.1716 10.000000 2
8 2 2 107.8284 10.000000 2
9 3 3 105.0000 7.171573 2
10 4 4 105.0000 12.828427 2
11 5 5 105.0000 10.000000 2
12 6 6 105.0000 10.000000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
#CCI
h.cci <- ccd(basis = 2, randomize = F, n0 = c(2,2), inscribed = T,
coding = list(x1~(la-105)/2,x2~(aa-10)/2),)
h.cci
run.order std.order la aa Block
1 1 1 103.5858 8.585786 1
2 2 2 106.4142 8.585786 1
3 3 3 103.5858 11.414214 1
4 4 4 106.4142 11.414214 1
5 5 5 105.0000 10.000000 1
6 6 6 105.0000 10.000000 1
7 1 1 103.0000 10.000000 2
8 2 2 107.0000 10.000000 2
9 3 3 105.0000 8.000000 2
10 4 4 105.0000 12.000000 2
11 5 5 105.0000 10.000000 2
12 6 6 105.0000 10.000000 2
Data are stored in coded form using these coding formulas ...
x1 ~ (la - 105)/2
x2 ~ (aa - 10)/2
# Escritura y lectura de un archivo
write.csv(x = h.ccf,file = "ccf.csv")
cccr <- read.csv(file = "ccf.csv")
library (readxl)
ejemplo <- read_excel("ejemplo.xlsx")
#install.packages("foreign")
library (foreign)
#read.mtp(file)