Superficie de respuesta I
Kenneth Roy Cabrera Torres
Martes 24 de abril de 2018
library(readxl)
library(rsm)
library(car)## Loading required package: carDatalibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericLectura de la base de datos
library(readxl)
rqca1 <- read_excel("reaccionQuimicaI.xlsx")
rqca1## # A tibble: 7 x 3
## Tiempo Temp Rendimiento
## <dbl> <dbl> <dbl>
## 1 80. 170. 80.5
## 2 80. 180. 81.5
## 3 90. 170. 82.0
## 4 90. 180. 83.5
## 5 85. 175. 83.9
## 6 85. 175. 84.3
## 7 85. 175. 84.0library(rsm)
# Codificación de las variables de acuerdo con la
# metodologÃa de superficie de respuesta.
rq1C <- coded.data(rqca1, x1 ~ (Tiempo-85)/5,
x2 ~ (Temp-175)/5)
rq1C## Tiempo Temp Rendimiento
## 1 80 170 80.5
## 2 80 180 81.5
## 3 90 170 82.0
## 4 90 180 83.5
## 5 85 175 83.9
## 6 85 175 84.3
## 7 85 175 84.0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Tiempo - 85)/5
## x2 ~ (Temp - 175)/5# Muestra la codificación correspondiente
as.data.frame(rq1C)## x1 x2 Rendimiento
## 1 -1 -1 80.5
## 2 -1 1 81.5
## 3 1 -1 82.0
## 4 1 1 83.5
## 5 0 0 83.9
## 6 0 0 84.3
## 7 0 0 84.0Modelización
# Modelo inicial de primer orden (FO = "Fisrt Order")
# Otros modelos son:
# TWI : (two-way interaction). Interacción doble.
# PQ : (pure cuadratic). Cuadrático puro.
# SO : (second-order). Segundo orden = FO + TWI + PQFO:
\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon \]
TWI:
\[ y = \beta_0 + \beta_3 x_1 x_2 + \epsilon \] PQ:
\[ y = \beta_0 + \beta_4 x_1^2 + \beta_5 x_2^2 + \epsilon \] SO:
\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 + \beta_4 x_1^2 + \beta_5 x_2^2 + \epsilon \]
m1.rq1 <- rsm(Rendimiento ~ FO(x1, x2), data = rq1C)
summary(m1.rq1)##
## Call:
## rsm(formula = Rendimiento ~ FO(x1, x2), data = rq1C)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.81429 0.54719 151.3456 1.143e-08 ***
## x1 0.87500 0.72386 1.2088 0.2933
## x2 0.62500 0.72386 0.8634 0.4366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.3555, Adjusted R-squared: 0.0333
## F-statistic: 1.103 on 2 and 4 DF, p-value: 0.4153
##
## Analysis of Variance Table
##
## Response: Rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 4.6250 2.3125 1.1033 0.41534
## Residuals 4 8.3836 2.0959
## Lack of fit 2 8.2969 4.1485 95.7335 0.01034
## Pure error 2 0.0867 0.0433
##
## Direction of steepest ascent (at radius 1):
## x1 x2
## 0.8137335 0.5812382
##
## Corresponding increment in original units:
## Tiempo Temp
## 4.068667 2.906191# Dirección de la mayor pendiente
steepest(m1.rq1)## Path of steepest ascent from ridge analysis:## dist x1 x2 | Tiempo Temp | yhat
## 1 0.0 0.000 0.000 | 85.000 175.000 | 82.814
## 2 0.5 0.407 0.291 | 87.035 176.455 | 83.352
## 3 1.0 0.814 0.581 | 89.070 177.905 | 83.890
## 4 1.5 1.221 0.872 | 91.105 179.360 | 84.428
## 5 2.0 1.627 1.162 | 93.135 180.810 | 84.964
## 6 2.5 2.034 1.453 | 95.170 182.265 | 85.502
## 7 3.0 2.442 1.744 | 97.210 183.720 | 86.041
## 8 3.5 2.848 2.035 | 99.240 185.175 | 86.578
## 9 4.0 3.255 2.325 | 101.275 186.625 | 87.116
## 10 4.5 3.662 2.616 | 103.310 188.080 | 87.654
## 11 5.0 4.069 2.906 | 105.345 189.530 | 88.191# Segundo modelo
m2.rq1 <- rsm(Rendimiento ~ FO(x1, x2) + TWI(x1, x2), data = rq1C)
summary(m2.rq1)##
## Call:
## rsm(formula = Rendimiento ~ FO(x1, x2) + TWI(x1, x2), data = rq1C)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.81429 0.62948 131.5604 9.683e-07 ***
## x1 0.87500 0.83272 1.0508 0.3705
## x2 0.62500 0.83272 0.7506 0.5074
## x1:x2 0.12500 0.83272 0.1501 0.8902
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.3603, Adjusted R-squared: -0.2793
## F-statistic: 0.5633 on 3 and 3 DF, p-value: 0.6755
##
## Analysis of Variance Table
##
## Response: Rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 4.6250 2.3125 0.8337 0.515302
## TWI(x1, x2) 1 0.0625 0.0625 0.0225 0.890202
## Residuals 3 8.3211 2.7737
## Lack of fit 1 8.2344 8.2344 190.0247 0.005221
## Pure error 2 0.0867 0.0433
##
## Stationary point of response surface:
## x1 x2
## -5 -7
##
## Stationary point in original units:
## Tiempo Temp
## 60 140
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.0625 -0.0625
##
## $vectors
## [,1] [,2]
## x1 0.7071068 -0.7071068
## x2 0.7071068 0.7071068# Gráfica para el segundo modelo
contour(m2.rq1, ~ x1 + x2, image = T, at = summary(m2.rq1)$canonical$xs)
# Dirección de la mayor pendiente
steepest(m2.rq1)## Path of steepest ascent from ridge analysis:## dist x1 x2 | Tiempo Temp | yhat
## 1 0.0 0.000 0.000 | 85.000 175.000 | 82.814
## 2 0.5 0.402 0.297 | 87.010 176.485 | 83.367
## 3 1.0 0.795 0.606 | 88.975 178.030 | 83.949
## 4 1.5 1.183 0.923 | 90.915 179.615 | 84.563
## 5 2.0 1.565 1.246 | 92.825 181.230 | 85.206
## 6 2.5 1.943 1.573 | 94.715 182.865 | 85.880
## 7 3.0 2.318 1.905 | 96.590 184.525 | 86.585
## 8 3.5 2.690 2.239 | 98.450 186.195 | 87.320
## 9 4.0 3.060 2.576 | 100.300 187.880 | 88.087
## 10 4.5 3.429 2.915 | 102.145 189.575 | 88.886
## 11 5.0 3.795 3.255 | 103.975 191.275 | 89.713Toma de datos de la segunda etapa (Segundo nivel del bloque)
rqca2 <- read_excel("reaccionQuimicaII.xlsx")
rqca2## # A tibble: 7 x 3
## Tiempo Temp Rendimiento
## <dbl> <dbl> <dbl>
## 1 85.0 175. 79.7
## 2 85.0 175. 79.8
## 3 85.0 175. 79.5
## 4 92.1 175. 78.4
## 5 77.9 175. 75.6
## 6 85.0 182. 78.5
## 7 85.0 168. 77.0# Unión de datos adicionales para mejorar el modelo
# mediante datos tomado en un segundo bloque
rq2C <- djoin(rq1C, rqca2)
rq2C## Tiempo Temp Rendimiento Block
## 1 80.00 170.00 80.5 1
## 2 80.00 180.00 81.5 1
## 3 90.00 170.00 82.0 1
## 4 90.00 180.00 83.5 1
## 5 85.00 175.00 83.9 1
## 6 85.00 175.00 84.3 1
## 7 85.00 175.00 84.0 1
## 8 85.00 175.00 79.7 2
## 9 85.00 175.00 79.8 2
## 10 85.00 175.00 79.5 2
## 11 92.07 175.00 78.4 2
## 12 77.93 175.00 75.6 2
## 13 85.00 182.07 78.5 2
## 14 85.00 167.93 77.0 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Tiempo - 85)/5
## x2 ~ (Temp - 175)/5as.data.frame(rq2C)## x1 x2 Rendimiento Block
## 1 -1.000 -1.000 80.5 1
## 2 -1.000 1.000 81.5 1
## 3 1.000 -1.000 82.0 1
## 4 1.000 1.000 83.5 1
## 5 0.000 0.000 83.9 1
## 6 0.000 0.000 84.3 1
## 7 0.000 0.000 84.0 1
## 8 0.000 0.000 79.7 2
## 9 0.000 0.000 79.8 2
## 10 0.000 0.000 79.5 2
## 11 1.414 0.000 78.4 2
## 12 -1.414 0.000 75.6 2
## 13 0.000 1.414 78.5 2
## 14 0.000 -1.414 77.0 2# Tercer modelo con los bloques y de segundo orden
m3.rq2 <- rsm(Rendimiento ~ Block + SO(x1, x2), data = rq2C)
summary(m3.rq2)##
## Call:
## rsm(formula = Rendimiento ~ Block + SO(x1, x2), data = rq2C)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.095427 0.079631 1056.067 < 2.2e-16 ***
## Block2 -4.457530 0.087226 -51.103 2.877e-10 ***
## x1 0.932541 0.057699 16.162 8.444e-07 ***
## x2 0.577712 0.057699 10.012 2.122e-05 ***
## x1:x2 0.125000 0.081592 1.532 0.1694
## x1^2 -1.308555 0.060064 -21.786 1.083e-07 ***
## x2^2 -0.933442 0.060064 -15.541 1.104e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9981, Adjusted R-squared: 0.9964
## F-statistic: 607.2 on 6 and 7 DF, p-value: 3.811e-09
##
## Analysis of Variance Table
##
## Response: Rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 1 69.531 69.531 2611.0950 2.879e-10
## FO(x1, x2) 2 9.626 4.813 180.7341 9.450e-07
## TWI(x1, x2) 1 0.063 0.063 2.3470 0.1694
## PQ(x1, x2) 2 17.791 8.896 334.0539 1.135e-07
## Residuals 7 0.186 0.027
## Lack of fit 3 0.053 0.018 0.5307 0.6851
## Pure error 4 0.133 0.033
##
## Stationary point of response surface:
## x1 x2
## 0.3722954 0.3343802
##
## Stationary point in original units:
## Tiempo Temp
## 86.86148 176.67190
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.9233027 -1.3186949
##
## $vectors
## [,1] [,2]
## x1 -0.1601375 -0.9870947
## x2 -0.9870947 0.1601375Gráfica del modelo final
contour(m3.rq2, ~ x1 + x2, image = T, at = summary(m3.rq2)$canonical$xs,
las = 1)
# Dirección de la mayor pendiente
steepest(m3.rq2)## Path of steepest ascent from ridge analysis:## dist x1 x2 | Tiempo Temp | yhat
## 1 0.0 0.000 0.000 | 85.000 175.000 | 84.095
## 2 0.5 0.372 0.334 | 86.860 176.670 | 84.366
## 3 1.0 0.640 0.768 | 88.200 178.840 | 84.111
## 4 1.5 0.838 1.244 | 89.190 181.220 | 83.362
## 5 2.0 0.995 1.735 | 89.975 183.675 | 82.136
## 6 2.5 1.129 2.230 | 90.645 186.150 | 80.441
## 7 3.0 1.249 2.728 | 91.245 188.640 | 78.274
## 8 3.5 1.360 3.225 | 91.800 191.125 | 75.646
## 9 4.0 1.465 3.723 | 92.325 193.615 | 72.548
## 10 4.5 1.564 4.219 | 92.820 196.095 | 69.000
## 11 5.0 1.660 4.716 | 93.300 198.580 | 64.980shapiro.test(residuals(m3.rq2))##
## Shapiro-Wilk normality test
##
## data: residuals(m3.rq2)
## W = 0.95814, p-value = 0.6924shapiro.test(rstandard(m3.rq2))##
## Shapiro-Wilk normality test
##
## data: rstandard(m3.rq2)
## W = 0.90573, p-value = 0.1366shapiro.test(rstudent(m3.rq2))##
## Shapiro-Wilk normality test
##
## data: rstudent(m3.rq2)
## W = 0.9245, p-value = 0.2551library(lmtest)
bptest(m3.rq2)##
## studentized Breusch-Pagan test
##
## data: m3.rq2
## BP = 7.759, df = 6, p-value = 0.2563library(lmtest)
bptest(m3.rq2, studentize = FALSE)##
## Breusch-Pagan test
##
## data: m3.rq2
## BP = 3.7093, df = 6, p-value = 0.7159library(car)
ncvTest(m3.rq2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.241426 Df = 1 p = 0.2651965op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(m3.rq2)
par(op)library(car)
influenceIndexPlot(m3.rq2)
predict(m3.rq2, interval = "confidence")## fit lwr upr
## 1 80.46818 80.14615 80.79021
## 2 81.37360 81.05157 81.69563
## 3 82.08326 81.76123 82.40529
## 4 83.48868 83.16665 83.81071
## 5 84.09543 83.90713 84.28372
## 6 84.09543 83.90713 84.28372
## 7 84.09543 83.90713 84.28372
## 8 79.63790 79.44962 79.82617
## 9 79.63790 79.44962 79.82617
## 10 79.63790 79.44962 79.82617
## 11 78.34019 78.01818 78.66219
## 12 75.70296 75.38096 76.02497
## 13 78.58846 78.26646 78.91047
## 14 76.95469 76.63269 77.27670predict(m3.rq2, interval = "prediction")## Warning in predict.lm(m3.rq2, interval = "prediction"): predictions on current data refer to _future_ responses## fit lwr upr
## 1 80.46818 79.96558 80.97077
## 2 81.37360 80.87101 81.87619
## 3 82.08326 81.58067 82.58585
## 4 83.48868 82.98609 83.99128
## 5 84.09543 83.66607 84.52479
## 6 84.09543 83.66607 84.52479
## 7 84.09543 83.66607 84.52479
## 8 79.63790 79.20855 80.06725
## 9 79.63790 79.20855 80.06725
## 10 79.63790 79.20855 80.06725
## 11 78.34019 77.83761 78.84277
## 12 75.70296 75.20039 76.20554
## 13 78.58846 78.08589 79.09104
## 14 76.95469 76.45212 77.45727pred1 <- data.frame(Tiempo = c(86.86, 86.86) ,
Temp = c(176.67, 176.67),
Block = c("1", "2"))
rsm_pred1 <- coded.data(pred1,
x1 ~ (Tiempo-85)/5,
x2 ~ (Temp-175)/5)prediccion1 <- predict(m3.rq2, newdata = rsm_pred1, interval = "confidence")
cbind(pred1, prediccion1)## Tiempo Temp Block fit lwr upr
## 1 86.86 176.67 1 84.36561 84.17809 84.55312
## 2 86.86 176.67 2 79.90808 79.72057 80.09558prediccion2 <- predict(m3.rq2, newdata = rsm_pred1, interval = "prediction")
cbind(pred1, prediccion2)## Tiempo Temp Block fit lwr upr
## 1 86.86 176.67 1 84.36561 83.93658 84.79463
## 2 86.86 176.67 2 79.90808 79.47906 80.33709