Paper 10 Proyecto Final
Angélica Reyes
10/23/2021
PROBLEMA PAPER 10
https://doi.org/10.1080/23311932.2020.1782097
Se utilizó un diseño experimental CCD para estudiar el efecto de la temperatura (X1) y la actividad del agua aw (X2) sobre el crecimiento por parte de la especie de Aspergillus flavus en un medio de cultivo fortificado con Chile en polvo.
library(knitr)
library(rsm)
library(tidyr)
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterlibrary(ggpubr)## Loading required package: ggplot2Cada factor se examinó en cinco niveles diferentes (-1,414, -1, 0, +1, +1,414) con cinco réplicas en el punto central. Los valores codificados y los valores reales de cada factor se presentan en la Tabla 1. Se realizaron un total de 13 series experimentales por triplicado.
table1<-read.csv("https://raw.githubusercontent.com/yvschemistry06/Experimental-Design/main/tabla%201.csv")
names(table1)=c("Factores","Simbolo","-1.414","-1","0","+1","+1.414")
kable(table1,"simple",caption = "Tabla 1. Factores y sus niveles codificados")| Factores | Simbolo | -1.414 | -1 | 0 | +1 | +1.414 |
|---|---|---|---|---|---|---|
| Temperatura | X1 | 22.930 | 25.0 | 30.000 | 35.00 | 37.070 |
| aw | X2 | 0.885 | 0.9 | 0.935 | 0.97 | 0.984 |
Con los datos de cada experimento (Tabla 2) determine:
- modelo matemático que mejor se ajusta,
- punto óptimo de mayor crecimiento,
- modelo de superficie de respuesta.
Haga una comparación de sus resultados con los obtenidos en el artículo científico.
df<-read.csv("https://raw.githubusercontent.com/yvschemistry06/Experimental-Design/main/Tabla%202.csv")
df## id temp aw Y
## 1 1 30.00 0.93500 7.32
## 2 2 30.00 0.93500 7.27
## 3 3 30.00 0.93500 7.10
## 4 4 35.00 0.97000 7.07
## 5 5 30.00 0.98449 8.30
## 6 6 22.93 0.93500 7.10
## 7 7 30.00 0.93500 7.23
## 8 8 30.00 0.93500 6.93
## 9 9 25.00 0.90000 4.23
## 10 10 35.00 0.90000 2.25
## 11 11 37.07 0.93500 3.97
## 12 12 25.00 0.97000 8.53
## 13 13 30.00 0.88551 3.10df <- coded.data(df, x1 ~ (temp - 30)/5, x2 ~ (aw - 0.935)/0.035)
str(df)## Classes 'coded.data' and 'data.frame': 13 obs. of 4 variables:
## $ id: int 1 2 3 4 5 6 7 8 9 10 ...
## $ x1: num 0 0 0 1 0 ...
## $ x2: num 0 0 0 1 1.41 ...
## $ Y : num 7.32 7.27 7.1 7.07 8.3 7.1 7.23 6.93 4.23 2.25 ...
## - attr(*, "codings")=List of 2
## ..$ x1:Class 'formula' language x1 ~ (temp - 30)/5
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## ..$ x2:Class 'formula' language x2 ~ (aw - 0.935)/0.035
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "rsdes")=List of 2
## ..$ primary: chr [1:2] "x1" "x2"
## ..$ call : language coded.data(data = df, x1 ~ (temp - 30)/5, x2 ~ (aw - 0.935)/0.035)Códigos
df$x1## [1] 0.000 0.000 0.000 1.000 0.000 -1.414 0.000 0.000 -1.000 1.000
## [11] 1.414 -1.000 0.000df$x2## [1] 0.000 0.000 0.000 1.000 1.414 0.000 0.000 0.000 -1.000 -1.000
## [11] 0.000 1.000 -1.414Ajuste de Primer Orden (Modelo lineal)
## df.rsm<-rsm(Yield~FO(x1,x2)+TWI(x1,x2)+PQ(x1,x2),data=df)
## df.rsm<-rsm(Yield~SO(x1,x2),data=df)
df.rsm.FO<-rsm(Y~FO(x1,x2),data=df)
summary(df.rsm.FO)##
## Call:
## rsm(formula = Y ~ FO(x1, x2), data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.18462 0.25740 24.0277 3.546e-10 ***
## x1 -0.98338 0.32814 -2.9968 0.01342 *
## x2 2.05941 0.32814 6.2760 9.190e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.8287, Adjusted R-squared: 0.7944
## F-statistic: 24.18 on 2 and 10 DF, p-value: 0.0001476
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 41.659 20.8297 24.184 0.0001476
## Residuals 10 8.613 0.8613
## Lack of fit 6 8.514 1.4190 57.567 0.0007744
## Pure error 4 0.099 0.0246
##
## Direction of steepest ascent (at radius 1):
## x1 x2
## -0.4308992 0.9024001
##
## Corresponding increment in original units:
## temp aw
## -2.154496 0.031584Lack of fit en el Pr debe ser Pr(>F)> 0.05, para que sea optima. el acontecido es bastante significativo con un p<0.0001 y un R-squared de 0.82. Para este modelo aún no obtenemos valores de predicción que satisfagan el diseño experimental.
Primer orden con interacción
df.INT<-rsm(Y~FO(x1,x2)+TWI(x1,x2),data=df)
summary(df.INT)##
## Call:
## rsm(formula = Y ~ FO(x1, x2) + TWI(x1, x2), data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.18462 0.27025 22.8847 2.759e-09 ***
## x1 -0.98338 0.34453 -2.8542 0.0189573 *
## x2 2.05941 0.34453 5.9774 0.0002082 ***
## x1:x2 0.13000 0.48720 0.2668 0.7956138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.83, Adjusted R-squared: 0.7734
## F-statistic: 14.65 on 3 and 9 DF, p-value: 0.000826
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 41.659 20.8297 21.9383 0.0003463
## TWI(x1, x2) 1 0.068 0.0676 0.0712 0.7956138
## Residuals 9 8.545 0.9495
## Lack of fit 5 8.447 1.6893 68.5322 0.0005758
## Pure error 4 0.099 0.0246
##
## Stationary point of response surface:
## x1 x2
## -15.841623 7.564431
##
## Stationary point in original units:
## temp aw
## -49.208114 1.199755
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.065 -0.065
##
## $vectors
## [,1] [,2]
## x1 0.7071068 -0.7071068
## x2 0.7071068 0.7071068Lack of fit del Pr(>F) 0.0005758 Para este modelo aún no obtenemos valores de predicción que satisfagan el diseño experimental.
Segundo orden sin interacción
df.rsm.PQ<-rsm(Y~FO(x1,x2)+PQ(x1,x2),data=df)
summary(df.rsm.PQ)##
## Call:
## rsm(formula = Y ~ FO(x1, x2) + PQ(x1, x2), data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.17001 0.13192 54.3498 1.457e-11 ***
## x1 -0.98338 0.10430 -9.4281 1.315e-05 ***
## x2 2.05941 0.10430 19.7446 4.507e-08 ***
## x1^2 -0.84202 0.11187 -7.5269 6.754e-05 ***
## x2^2 -0.75949 0.11187 -6.7892 0.0001393 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9862, Adjusted R-squared: 0.9792
## F-statistic: 142.4 on 4 and 8 DF, p-value: 1.818e-07
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 41.659 20.8297 239.3695 7.298e-08
## PQ(x1, x2) 2 7.917 3.9583 45.4881 4.268e-05
## Residuals 8 0.696 0.0870
## Lack of fit 4 0.598 0.1494 6.0604 0.0545
## Pure error 4 0.099 0.0246
##
## Stationary point of response surface:
## x1 x2
## -0.5839401 1.3557801
##
## Stationary point in original units:
## temp aw
## 27.0802996 0.9824523
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.7594930 -0.8420179
##
## $vectors
## [,1] [,2]
## x1 0 -1
## x2 -1 0Multiple R-squared: 0.9862 aumentó Lack of fit de Pr(>F) 0.0545 es no significativo Los valores partida codificados para realizar el experimento son: Stationary point of response surface: x1 x2 -0.5839401 1.3557801
R= Este modelo es el más óptimo, ya que su predicción se acerca a 1 con el R-squared obteniendo los valores reales para proyectar el experimento con temperatura 27.08 y Actividad del agua (aw) 0.9824523. Se debe realizar el ANOVA para comprobar las diferencias entre los modelos más Óptimos.
Modelo completo (Segundo orden)
df.rsm.SO<-rsm(Y~SO(x1,x2),data=df)
summary(df.rsm.SO)##
## Call:
## rsm(formula = Y ~ SO(x1, x2), data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.17001 0.13401 53.5037 2.088e-10 ***
## x1 -0.98338 0.10595 -9.2813 3.490e-05 ***
## x2 2.05941 0.10595 19.4372 2.380e-07 ***
## x1:x2 0.13000 0.14983 0.8677 0.4143302
## x1^2 -0.84202 0.11364 -7.4097 0.0001482 ***
## x2^2 -0.75949 0.11364 -6.6835 0.0002817 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9875, Adjusted R-squared: 0.9786
## F-statistic: 110.6 on 5 and 7 DF, p-value: 1.671e-06
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 41.659 20.8297 231.9743 4.003e-07
## TWI(x1, x2) 1 0.068 0.0676 0.7528 0.4143302
## PQ(x1, x2) 2 7.917 3.9583 44.0828 0.0001079
## Residuals 7 0.629 0.0898
## Lack of fit 3 0.530 0.1767 7.1663 0.0436512
## Pure error 4 0.099 0.0246
##
## Stationary point of response surface:
## x1 x2
## -0.4824674 1.3144889
##
## Stationary point in original units:
## temp aw
## 27.5876628 0.9810071
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.7237646 -0.8777463
##
## $vectors
## [,1] [,2]
## x1 -0.4816950 -0.8763389
## x2 -0.8763389 0.4816950En el intercepto x1:x2 Pr(>|t|) es no significativo ya que es mayor de alpha al 95% de confidencialidad
Multiple R-squared: 0.9786
Lack of fit 0.0436512 significativo, no supera el presentado en el Segundo orden sin interacción.
Se debe realizar el ANOVA para comprobar las diferencias entre los modelos más Óptimos.
ANOVA para ver si el aumento en la complejidad del modelo es significativa
anova(df.rsm.PQ,df.rsm.SO)## Analysis of Variance Table
##
## Model 1: Y ~ FO(x1, x2) + PQ(x1, x2)
## Model 2: Y ~ FO(x1, x2) + TWI(x1, x2) + PQ(x1, x2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 8 0.69615
## 2 7 0.62855 1 0.0676 0.7528 0.4143Es no significativo por lo que no es necesario incluir la interacción.
Intervalos de confianza
confint(df.rsm.PQ)## 2.5 % 97.5 %
## (Intercept) 6.865796 7.4742275
## FO(x1, x2)x1 -1.223898 -0.7428541
## FO(x1, x2)x2 1.818889 2.2999328
## PQ(x1, x2)x1^2 -1.099985 -0.5840505
## PQ(x1, x2)x2^2 -1.017460 -0.5015256#Predicción a partir del Modelo seleccionado
predict(df.rsm.PQ,df[1:12,],interval="confidence")## fit lwr upr
## 1 7.170012 6.865796 7.474227
## 2 7.170012 6.865796 7.474227
## 3 7.170012 6.865796 7.474227
## 4 6.644536 6.227929 7.061142
## 5 8.563492 8.025741 9.101242
## 6 6.876978 6.339228 7.414728
## 7 7.170012 6.865796 7.474227
## 8 7.170012 6.865796 7.474227
## 9 4.492466 4.075859 4.909072
## 10 2.525714 2.109107 2.942320
## 11 4.095991 3.558241 4.633741
## 12 8.611288 8.194681 9.027894Los valores reales para proyectar el experimento son el modelo de Segundo orden sin interacción (df.rsm.PQ) con temperatura 27.0802996 y Actividad del agua (aw) 0.9824523.
par(mfrow = c(1, 1))
contour(df.rsm.PQ, ~ x1 + x2, image = TRUE,xlabs=c("Actividad del Agua(aw)","Temperatura (°C)"))
abline(v=27.0802996,lty=3)
abline(h=0.9824523,lty=3)
text(27.0,0.979,"(27.0802996,0.9824523)")
persp(df.rsm.PQ,x1~x2,col=terrain.colors(50),
contours="color",
zlab="Y(%)",
xlabs=c("Actividad del Agua(aw)","Temperatura (°C)",
theta= 80, phi= 90)
)
Al observar la gráfica de superficie se nota que no tiene un diseño visual de terminado para los modelos de segundo orden como los gráficos regulares observados en el Capitulo 12 pág 350 de Gutierrez y de la Vara, 3ra ed. Análisis y Diseño experimentos. Sin embargo, el gráfico si muestra los datos introducidos y el punto predictivo como optimo según cálculos estadísticos. Según literatura, se puede probar con un análisis canónico o un análisis de cordillera dependiendo de lo que se desee.