Solicitud de la consultante
Conocer las mejores condiciones experimentales para obtener una pectina con un alto contenido de ácido galacturónico y grado de esterificación.
Metodología propuesta
- Diseño experimental con tres factores.
- Análisis de varianza multivariado (MANOVA)
- Metodología de Superficies de Respuesta (RSM)
Uso de software libre
Instalación de librerías
Se requiere de la instalación de las siguientes librerías, solamente debe ejecutarse este comando en la consola de R una sola vez y se instalarán las librerías respectivas:
install.packages("readxl")
install.packages("rsm")
install.packages("glmtoolbox")
install.packages("ggplot2")
install.packages("MVN")
install.packages("car")
install.packages("biotools")
install.packages("dplyr")
install.packages("knitr")Así mismo siempre debe hacerse el siguiente paso para llamar las librerías respectivas:
library(readxl)
library(rsm)
library(glmtoolbox)
library(ggplot2)
library(MVN)## Warning: package 'MVN' was built under R version 4.2.2library(car)## Loading required package: carDatalibrary(biotools)## Warning: package 'biotools' was built under R version 4.2.2## Loading required package: MASS## ---
## biotools version 4.2library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)Lectura de los datos del estudio
Para realizar el análisis debemos leer el conjunto de datos que usaremos en el estudio, en este caso se encuentra en MUCILAGO.xlsx, se envía el archivo como adjunto para replicar los resultados obtenidos en este script.
## Lectura de los datos
MUCILAGO <- read_excel("MUCILAGO.xlsx")
## Creación de inv_LS porque LS no está igualmente espaciada
MUCILAGO$inv_LS<-1/MUCILAGO$LS
knitr::kable(x = MUCILAGO, digits = 2, align = "c")| Ensayo | pH | t_min | LS | DE | GA | Yield | Yield_GA | inv_LS |
|---|---|---|---|---|---|---|---|---|
| 1 C | 2.5 | 7.5 | 0.07 | 78.61 | 40.85 | 10.81 | 4.42 | 15 |
| 2 | 2.5 | 5.0 | 0.10 | 65.36 | 48.48 | 8.13 | 3.94 | 10 |
| 3 | 3.0 | 10.0 | 0.07 | 64.28 | 36.81 | 8.70 | 3.20 | 15 |
| 4 | 2.0 | 7.5 | 0.05 | 80.18 | 61.07 | 14.76 | 9.01 | 20 |
| 5 C | 2.5 | 7.5 | 0.07 | 82.47 | 42.86 | 9.36 | 4.01 | 15 |
| 6 | 3.0 | 7.5 | 0.10 | 64.61 | 36.31 | 8.89 | 3.23 | 10 |
| 7 | 2.5 | 10.0 | 0.10 | 75.65 | 48.64 | 12.74 | 6.20 | 10 |
| 8 | 3.0 | 5.0 | 0.07 | 74.95 | 29.08 | 10.67 | 3.10 | 15 |
| 9 | 2.5 | 5.0 | 0.05 | 62.66 | 49.62 | 9.48 | 4.70 | 20 |
| 10 | 3.0 | 7.5 | 0.05 | 65.29 | 34.75 | 8.50 | 2.95 | 20 |
| 11 C | 2.5 | 7.5 | 0.07 | 77.25 | 48.90 | 9.98 | 4.88 | 15 |
| 12 | 2.5 | 10.0 | 0.05 | 69.19 | 44.98 | 11.68 | 5.26 | 20 |
| 13 | 2.0 | 10.0 | 0.07 | 77.14 | 61.89 | 12.49 | 7.73 | 15 |
| 14 C | 2.5 | 7.5 | 0.07 | 82.76 | 43.54 | 10.11 | 4.40 | 15 |
| 15 | 2.0 | 7.5 | 0.10 | 81.79 | 54.62 | 13.90 | 7.59 | 10 |
| 16 | 2.0 | 5.0 | 0.07 | 72.89 | 56.29 | 11.92 | 6.71 | 15 |
Nuestras variables de interés son el grado de esterificación, el porcentaje de ácido galacturónico y el rendimiento.
Análisis Descriptivo de los datos
Variables por separado
Analizamos el comportamiento de las variables por separado, comenzamos con el grado de esterificación a la luz de los tres factores explicativos
plot(as.factor(MUCILAGO$pH),MUCILAGO$DE, xlab="Niveles de pH", main= "pH vs. % Grado de Esterificación") ### pH vs. DE
plot(as.factor(MUCILAGO$t_min),MUCILAGO$DE, xlab="Niveles del tiempo", main= "Tiempo de extracción vs. % Grado de Esterificación") ### tiempo vs. DE
plot(as.factor(MUCILAGO$inv_LS),MUCILAGO$DE, xlab="Niveles de v/m",main= "v/m vs. % Grado de Esterificación") ### v/m vs. DE
Estos gráficos nos indican cómo se comportan el grado de esterificación
respecto a las tres variables explicativas y nos generan indicios de lo
que puede presentarse con cada una por separado. En el caso del pH se
ven diferencias entre las medianas de cada nivel de pH, razón que la
hace candidata a ser significativa en el modelo, sin embargo la
dispersión presente en las otras dos variables (t_min y v/m) es un
indicio de que tal vez no aporten mucha información al modelo del
rendimiento
Continuamos con el porcentaje de ácido galacturónico explicado por los tres factores
plot(as.factor(MUCILAGO$pH),MUCILAGO$GA,main= "pH vs. % Ácido Galacturónico") ### pH vs. GA
plot(as.factor(MUCILAGO$t_min),MUCILAGO$GA,main= "Tiempo de extracción vs. % Ácido Galacturónico") ### tiempo vs. GA
plot(as.factor(MUCILAGO$inv_LS),MUCILAGO$GA,main= "Relación v/m vs. % Ácido Galacturónico") ### v/m vs. GA
Respecto al ácido galacturónico es evidente que hay diferencias
significativas entre los niveles de pH, mientras que para las otras
variables se ve un comportamiento muy parecido entre las medianas de
cada nivel de las variables.
Analizamos la última variable del estudio
plot(as.factor(MUCILAGO$pH),MUCILAGO$Yield_GA, xlab="Niveles de pH",main= "pH vs. Rendimiento (GA)") ### pH vs. Rendimiento
plot(as.factor(MUCILAGO$t_min),MUCILAGO$Yield_GA, xlab="Niveles de tiempo",main= "Tiempo de Extracción vs. Rendimiento (GA)") ### tiempos vs. Rendimiento
plot(as.factor(MUCILAGO$inv_LS),MUCILAGO$Yield_GA, xlab="Niveles de m/v",main= "Relación v/m vs. Rendimiento (GA)") ### v/m vs. Rendimiento
Al igual que con las otras dos variables explicativas notamos que el pH
presenta indicios de tener diferencias en el rendimiento dependiendo de
su nivel específico, mientras que esto no es evidente con el tiempo de
extracción y la relación masa/volumen
Variables 2-2
DE explicado por pH vs inv_LS
ggplot(MUCILAGO,aes(x=factor(pH),y=DE,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Grado de esterificación explicado por relación v/m y el pH",x="Relación v/m", y="Grado de esterificación", fill= "Nivel de \n pH") +
facet_wrap(vars(inv_LS),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
De este gráfico podemos ver que quizá haya una interacción en el efecto de la relación m/v y el efecto del pH sobre el grado de esterificación
DE explicado por pH vs t_min
ggplot(MUCILAGO,aes(x=factor(pH),y=DE,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Grado de esterificación explicado por tiempo de extracción y el pH",x="Tiempo de extracción", y="Grado de esterificación", fill= "Nivel de \n pH") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
En este gráfico encontramos que es importante observar la relación existente entre los niveles del tiempo de extracción respecto a los niveles de pH, puesto que pareciera que tienen un efecto sobre el grado de esterificación.
DE explicado por inv_LS vs t_min
ggplot(MUCILAGO,aes(x=factor(inv_LS),y=DE,fill=factor(inv_LS))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Grado de esterificación explicado por tiempo de extracción y relación v/m",x="Tiempo de extracción", y="Grado de esterificación", fill= "Nivel de \n Relación m/v") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
De aquí observamos que tal vez el tiempo de extracción y el nivel de la relación de m/v tengan relación respecto al efecto conjunto que tiempo sobre el grado de esterificación.
GA explicado por pH vs inv_LS
ggplot(MUCILAGO,aes(x=factor(pH),y=GA,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="% ácido galacturónico explicado por relación mv/v y el pH",x="Relación m/v", y="% Ácido Galacturónico", fill= "Nivel de \n pH") +
facet_wrap(vars(inv_LS),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
Para este gráfico no se ven diferencias significativas en el efecto del pH al cambiar de nivel de la relación m/v
GA explicado por pH vs t_min
ggplot(MUCILAGO,aes(x=factor(pH),y=GA,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="% ácido galacturónico explicado por tiempo de extracción y el pH",x="Tiempo de extracción", y="% Ácido Galacturónico", fill= "Nivel de \n pH") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
En este gráfico tampoco se ven diferencias en el efecto del pH sobre el ácido galacturónico al cambiar el tiempo de extracción.
GA explicado por inv_LS vs t_min
ggplot(MUCILAGO,aes(x=factor(inv_LS),y=GA,fill=factor(inv_LS))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="% ácido galacturónico explicado por tiempo de extracción y la relación v/m",x="Tiempo de extracción", y="% Ácido Galacturónico", fill= "Nivel de \n Relación v/m") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
Acá tampoco hay diferencias significativas entre el efecto de la relación v/m sobre el ácido galacturónico al cambiar de nivel del tiempo de extracción.
Yield_GA explicado por pH vs inv_LS
ggplot(MUCILAGO,aes(x=factor(pH),y=Yield_GA,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Rendimiento explicado por pH y la relación v/m",x="Relación v/m", y="Rendimiento", fill= "Nivel de \n pH") +
facet_wrap(vars(inv_LS),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
No se observan diferencias en los efectos del pH al cambiar de nivel en la relación m/v
Yield_GA explicado por pH vs t_min
ggplot(MUCILAGO,aes(x=factor(pH),y=Yield_GA,fill=factor(pH))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Rendimiento explicado por pH y el tiempo de extracción",x="Tiempo de extracción", y="Rendimiento", fill= "Nivel de \n pH") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
No hay diferencias en el rendimiento respecto al efecto del pH cuando cambia el tiempo de extracción.
Yield_GA explicado por inv_LS vs t_min
ggplot(MUCILAGO,aes(x=factor(inv_LS),y=Yield_GA,fill=factor(inv_LS))) +
geom_boxplot(outlier.shape=NA,color="black",linetype="solid") +
scale_fill_manual(values=c("red","yellow","blue")) +
labs(title="Rendimiento explicado por el tiempo de extracción y la relación v/m",x="Tiempo de extracción", y="Rendimiento", fill= "Nivel de \n Relación v/m") +
facet_wrap(vars(t_min),ncol=3,nrow=1,strip.position="top",dir="h",scales="fixed")
En este gráfico tampoco se denotan diferencias en el rendimiento
respecto al efecto de la relación v/m al variar el tiempo de
extracción.
Análisis de varianza multivariado
Como es de interés conocer los valores que optimizan en conjunto a las variables DE y GA, en esta sección realizaremos un Análisis Multivariado para ejecutar un modelo multivariado y verificar si hay alguna variable que explique el comportamiento conjunto de las variables respuesta.
Verificamos supuestos
**Normalidad Multivariada y Normalidad Univariada+*
Normal_MUCILAGO<-mvn(MUCILAGO[,c(5,6)],mvnTest = "hz")
Normal_MUCILAGO$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.3992475 0.3971698 YESNormal_MUCILAGO$univariateNormality## Test Variable Statistic p value Normality
## 1 Anderson-Darling DE 0.6014 0.0980 YES
## 2 Anderson-Darling GA 0.1741 0.9097 YESLas variables en conjunto siguen una distribución normal multivariada y por aparte también se distribuyen normal.
Homogeneidad de varianzas multivariada
Presentamos primero la homogeneidad de varianzas para las 3 variables explicativas.
boxM(MUCILAGO[,c(5,6)],grouping = MUCILAGO$pH) ## P. de Homogeneidad respecto al pH##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: MUCILAGO[, c(5, 6)]
## Chi-Sq (approx.) = 7.8744, df = 6, p-value = 0.2475boxM(MUCILAGO[,c(5,6)],grouping = MUCILAGO$t_min) ## P. de Homogeneidad respecto al T. de extracción##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: MUCILAGO[, c(5, 6)]
## Chi-Sq (approx.) = 5.5007, df = 6, p-value = 0.4814boxM(MUCILAGO[,c(5,6)],grouping = MUCILAGO$inv_LS) ## P. de Homogeneidad respecto a la relación v/m##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: MUCILAGO[, c(5, 6)]
## Chi-Sq (approx.) = 2.5838, df = 6, p-value = 0.859No hay razón para creer que hay diferencias de varianzas entre los grupos elegidos.
Homogeneidad de varianzas univariada
Variable rendimiento en unidades de GA:
bartlett.test(MUCILAGO$Yield_GA~MUCILAGO$pH)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$Yield_GA by MUCILAGO$pH
## Bartlett's K-squared = 7.2638, df = 2, p-value = 0.02647bartlett.test(MUCILAGO$Yield_GA~MUCILAGO$t_min)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$Yield_GA by MUCILAGO$t_min
## Bartlett's K-squared = 0.36349, df = 2, p-value = 0.8338bartlett.test(MUCILAGO$Yield_GA~MUCILAGO$inv_LS)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$Yield_GA by MUCILAGO$inv_LS
## Bartlett's K-squared = 0.82701, df = 2, p-value = 0.6613Variable porcentaje de ácido galacturónico:
bartlett.test(MUCILAGO$GA~MUCILAGO$pH)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$GA by MUCILAGO$pH
## Bartlett's K-squared = 0.022224, df = 2, p-value = 0.9889bartlett.test(MUCILAGO$GA~MUCILAGO$t_min)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$GA by MUCILAGO$t_min
## Bartlett's K-squared = 0.2833, df = 2, p-value = 0.8679bartlett.test(MUCILAGO$GA~MUCILAGO$inv_LS)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$GA by MUCILAGO$inv_LS
## Bartlett's K-squared = 0.38741, df = 2, p-value = 0.8239Variable grado de esterificación:
bartlett.test(MUCILAGO$DE~MUCILAGO$pH)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$DE by MUCILAGO$pH
## Bartlett's K-squared = 1.6268, df = 2, p-value = 0.4433bartlett.test(MUCILAGO$DE~MUCILAGO$t_min)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$DE by MUCILAGO$t_min
## Bartlett's K-squared = 0.29675, df = 2, p-value = 0.8621bartlett.test(MUCILAGO$DE~MUCILAGO$inv_LS)##
## Bartlett test of homogeneity of variances
##
## data: MUCILAGO$DE by MUCILAGO$inv_LS
## Bartlett's K-squared = 0.55275, df = 2, p-value = 0.7585Ejecución del MANOVA
Modelo saturado
Se toma el modelo saturado para observar cuáles son los efectos más significativos.
modelo_conjunto<-manova(cbind(DE,GA)~inv_LS*pH*t_min,data=MUCILAGO)
summary(modelo_conjunto)## Df Pillai approx F num Df den Df Pr(>F)
## inv_LS 1 0.03154 0.130 2 8 0.8797
## pH 1 0.96705 117.391 2 8 1.179e-06 ***
## t_min 1 0.28470 1.592 2 8 0.2618
## inv_LS:pH 1 0.18865 0.930 2 8 0.4334
## inv_LS:t_min 1 0.15434 0.730 2 8 0.5114
## pH:t_min 1 0.15047 0.709 2 8 0.5208
## Residuals 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo 2
Se toman los elementos que aún se consideran significativos y se vuelve a evaluar
m2<-manova(cbind(DE,GA)~pH*t_min,data=MUCILAGO)
summary(m2)## Df Pillai approx F num Df den Df Pr(>F)
## pH 1 0.95422 114.639 2 11 4.303e-08 ***
## t_min 1 0.22289 1.577 2 11 0.2498
## pH:t_min 1 0.13348 0.847 2 11 0.4548
## Residuals 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo 3
Se siguen eliminando los efectos no significativos.
m3<-manova(cbind(DE,GA)~pH+t_min,data=MUCILAGO)
summary(m3)## Df Pillai approx F num Df den Df Pr(>F)
## pH 1 0.95189 118.720 2 12 1.24e-08 ***
## t_min 1 0.21005 1.595 2 12 0.243
## Residuals 13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo Final
Nos quedamos con el modelo univariado que explica la relación solo por medio del pH
mF<-manova(cbind(DE,GA)~pH,data=MUCILAGO)
summary(mF)## Df Pillai approx F num Df den Df Pr(>F)
## pH 1 0.94028 102.34 2 13 1.109e-08 ***
## Residuals 14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Superficies de respuesta univariadas
Codificación para las superficies de Respuesta
MUCILAGO1 <- coded.data(
MUCILAGO,
formulas = list(
x1 ~ (pH - 2.5)/0.5,
x2 ~ (t_min - 7.5)/2.5,
x3 ~ (inv_LS - 15)/5
))Modelo para el Rendimiento (Yield_GA)
mc1<- rsm(Yield_GA ~ SO(x1, x2, x3), data = MUCILAGO1)
summary(mc1)##
## Call:
## rsm(formula = Yield_GA ~ SO(x1, x2, x3), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.427364 0.258687 17.1147 2.547e-06 ***
## x1 -2.320167 0.182920 -12.6841 1.472e-05 ***
## x2 0.490892 0.182920 2.6836 0.03636 *
## x3 0.120937 0.182920 0.6611 0.53307
## x1:x2 -0.230298 0.258687 -0.8903 0.40760
## x1:x3 -0.423791 0.258687 -1.6382 0.15249
## x2:x3 -0.426016 0.258687 -1.6468 0.15069
## x1^2 0.715772 0.258687 2.7669 0.03255 *
## x2^2 0.043675 0.258687 0.1688 0.87148
## x3^2 0.553404 0.258687 2.1393 0.07623 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9689, Adjusted R-squared: 0.9223
## F-statistic: 20.77 on 9 and 6 DF, p-value: 0.0007422
##
## Analysis of Variance Table
##
## Response: Yield_GA
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 45.110 15.0367 56.1750 8.773e-05
## TWI(x1, x2, x3) 3 1.657 0.5522 2.0628 0.2067
## PQ(x1, x2, x3) 3 3.282 1.0940 4.0870 0.0673
## Residuals 6 1.606 0.2677
## Lack of fit 3 1.230 0.4100 3.2693 0.1783
## Pure error 3 0.376 0.1254
##
## Stationary point of response surface:
## x1 x2 x3
## 1.5302337 -0.8427949 0.1522559
##
## Stationary point in original units:
## pH t_min inv_LS
## 3.265117 5.393013 15.761279
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.8624687 0.5270875 -0.0767044
##
## $vectors
## [,1] [,2] [,3]
## x1 0.81201270 0.5356169 -0.2318403
## x2 0.03732681 -0.4440780 -0.8952103
## x3 -0.58244491 0.7182683 -0.3805898Veamos la gráfica de la superficie de respuesta del modelo propuesto inicialmente:
par(mfrow=c(1,3))
contour (mc1, ~ x1 + x2 + x3, image=TRUE)
persp(mc1,~x1+x2+x3, col="lightblue",contours = list(z="top", col="orange"),theta = -145, phi = 35)
Eliminamos los efectos que generan ruido sobre el modelo y obtenemos el siguiente modelo
mc1.1<- rsm(Yield_GA ~ FO(x1, x2)+PQ(x1), data = MUCILAGO1)
summary(mc1.1)## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend##
## Call:
## rsm(formula = Yield_GA ~ FO(x1, x2) + PQ(x1), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.72590 0.21919 21.5608 5.791e-11 ***
## x1 -2.32017 0.21919 -10.5852 1.931e-07 ***
## x2 0.49089 0.21919 2.2396 0.04483 *
## x1^2 0.71577 0.30998 2.3091 0.03954 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9107, Adjusted R-squared: 0.8884
## F-statistic: 40.8 on 3 and 12 DF, p-value: 1.428e-06
##
## Analysis of Variance Table
##
## Response: Yield_GA
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 44.993 22.4966 58.5312 6.461e-07
## PQ(x1) 1 2.049 2.0493 5.3319 0.03954
## Residuals 12 4.612 0.3844
## Lack of fit 5 2.456 0.4912 1.5945 0.27699
## Pure error 7 2.156 0.3080
##
## Stationary point of response surface:
## x1 x2
## 1.620744 0.000000
##
## Stationary point in original units:
## pH t_min
## 3.310372 7.500000
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.7157723 0.0000000
##
## $vectors
## [,1] [,2]
## x1 -1 0
## x2 0 -1Debemos verificar los supuestos del modelo para poder usarlo
shapiro.test(resid(mc1.1)) ## Normalidad##
## Shapiro-Wilk normality test
##
## data: resid(mc1.1)
## W = 0.91258, p-value = 0.1281bartlett.test(resid(mc1.1)~MUCILAGO1$x1) ### Homoscedasticidad##
## Bartlett test of homogeneity of variances
##
## data: resid(mc1.1) by MUCILAGO1$x1
## Bartlett's K-squared = 1.8541, df = 2, p-value = 0.3957bartlett.test(resid(mc1.1)~MUCILAGO1$x2)##
## Bartlett test of homogeneity of variances
##
## data: resid(mc1.1) by MUCILAGO1$x2
## Bartlett's K-squared = 0.18344, df = 2, p-value = 0.9124Se cumplen los supuestos y el modelo es acorde al ajuste. De esta manera obtenemos la siguiente superficie de nivel que explica todos los efectos significativos en el modelo y permite identificar el punto donde se optimiza el rendimiento.
contour (mc1.1, ~ x1 + x2, image=TRUE)
persp(mc1.1,~x1+x2, col="lightblue",contours = list(z="top", col="orange"),theta = -210, phi = 35)
Modelo para el Grado de Esterificación (DE)
Generamos un modelo inicial
mc2<- rsm(DE ~ SO(x1, x2, x3), data = MUCILAGO1)
summary(mc2)##
## Call:
## rsm(formula = DE ~ SO(x1, x2, x3), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.2725 2.5450 31.5415 6.748e-08 ***
## x1 -5.3587 1.7996 -2.9778 0.02471 *
## x2 1.3000 1.7996 0.7224 0.49723
## x3 -1.2612 1.7996 -0.7009 0.50964
## x1:x2 -3.7300 2.5450 -1.4656 0.19310
## x1:x3 0.5725 2.5450 0.2250 0.82948
## x2:x3 -0.9400 2.5450 -0.3694 0.72455
## x1^2 -1.6025 2.5450 -0.6297 0.55212
## x2^2 -6.3550 2.5450 -2.4971 0.04671 *
## x3^2 -5.7025 2.5450 -2.2407 0.06628 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.7991, Adjusted R-squared: 0.4978
## F-statistic: 2.652 on 9 and 6 DF, p-value: 0.1239
##
## Analysis of Variance Table
##
## Response: DE
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 255.976 85.325 3.2934 0.09975
## TWI(x1, x2, x3) 3 60.497 20.166 0.7784 0.54754
## PQ(x1, x2, x3) 3 301.890 100.630 3.8842 0.07408
## Residuals 6 155.446 25.908
## Lack of fit 3 132.530 44.177 5.7833 0.09172
## Pure error 3 22.916 7.639
##
## Stationary point of response surface:
## x1 x2 x3
## -2.8540466 0.9645115 -0.3333479
##
## Stationary point in original units:
## pH t_min inv_LS
## 1.072977 9.911279 13.333261
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.9202607 -5.6496937 -7.0900456
##
## $vectors
## [,1] [,2] [,3]
## x1 0.93974919 -0.1641962 -0.2998518
## x2 -0.33015760 -0.2083314 -0.9206487
## x3 0.08869846 0.9641772 -0.2499898Veamos la gráfica de la superficie de respuesta del modelo inicial, con variables de ruido:
par(mfrow=c(1,3))
contour (mc2, ~ x1 + x2 + x3, image=TRUE)
persp(mc2,~x1+x2+x3, col="lightblue",contours = list(z="top", col="orange"),theta = -145, phi = 35)
Eliminamos los efectos que generan ruido sobre el modelo y obtenemos el siguiente modelo
mc2.1<- rsm(DE ~ FO(x1, x2)+PQ(x2), data = MUCILAGO1)
summary(mc2.1)## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend##
## Call:
## rsm(formula = DE ~ FO(x1, x2) + PQ(x2), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.6200 1.9606 39.0801 5.082e-14 ***
## x1 -5.3588 1.9606 -2.7332 0.01816 *
## x2 1.3000 1.9606 0.6631 0.51982
## x2^2 -6.3550 2.7727 -2.2920 0.04078 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.5231, Adjusted R-squared: 0.4039
## F-statistic: 4.388 on 3 and 12 DF, p-value: 0.02649
##
## Analysis of Variance Table
##
## Response: DE
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 243.25 121.625 3.9551 0.047933
## PQ(x2) 1 161.54 161.544 5.2532 0.040784
## Residuals 12 369.02 30.751
## Lack of fit 5 320.06 64.012 9.1532 0.005606
## Pure error 7 48.95 6.993
##
## Stationary point of response surface:
## x1 x2
## 0.0000000 0.1022817
##
## Stationary point in original units:
## pH t_min
## 2.500000 7.755704
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.000 -6.355
##
## $vectors
## [,1] [,2]
## x1 -1 0
## x2 0 1Verificamos que el modelo corresponda al ajustado, verificando supuestos.
shapiro.test(resid(mc2.1))##
## Shapiro-Wilk normality test
##
## data: resid(mc2.1)
## W = 0.94713, p-value = 0.4457bartlett.test(resid(mc2.1)~MUCILAGO1$x1)##
## Bartlett test of homogeneity of variances
##
## data: resid(mc2.1) by MUCILAGO1$x1
## Bartlett's K-squared = 8.3012, df = 2, p-value = 0.01576bartlett.test(resid(mc2.1)~(MUCILAGO1$x2^2))##
## Bartlett test of homogeneity of variances
##
## data: resid(mc2.1) by MUCILAGO1$x2
## Bartlett's K-squared = 2.4839, df = 2, p-value = 0.2888Este modelo no es adecuado porque no cumple con uno de los supuestos. Ajustamos un modelo adicional para solucionar los problemas asociados con la varianza del modelo y encontramos que el Modelo Gamma explica muy bien la varianza de los datos y genera que se comporte bien respecto a sus residuales. Usando como función de enlace la función identidad tenemos que el análisis es análogo, pero en esta caso, podemos decir estadísticamente que la variable que explica el comportamiento del Grado de Esterificación es el pH.
fit_final<-glm(DE~pH, family = Gamma(identity),data=MUCILAGO) ##Modelo Gamma
summary(fit_final) ### Resumen del modelo##
## Call:
## glm(formula = DE ~ pH, family = Gamma(identity), data = MUCILAGO)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.154715 -0.056215 -0.003456 0.055393 0.121836
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.834 11.334 8.896 3.9e-07 ***
## pH -10.956 4.424 -2.476 0.0267 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.007297954)
##
## Null deviance: 0.14646 on 15 degrees of freedom
## Residual deviance: 0.10280 on 14 degrees of freedom
## AIC: 108
##
## Number of Fisher Scoring iterations: 4plot(resid(fit_final)) ### Gráfico de residuales, no hay tendencia
envelope(fit_final) ### Ajuste de los datos, se mantienen entre las bandas##
|
| | 0%
|
|+ | 4%
|
|++ | 8%
|
|+++ | 12%
|
|++++ | 16%
|
|+++++ | 20%
|
|++++++ | 24%
|
|+++++++ | 28%
|
|++++++++ | 32%
|
|+++++++++ | 36%
|
|++++++++++ | 40%
|
|+++++++++++ | 44%
|
|++++++++++++ | 48%
|
|+++++++++++++ | 52%
|
|++++++++++++++ | 56%
|
|+++++++++++++++ | 60%
|
|++++++++++++++++ | 64%
|
|+++++++++++++++++ | 68%
|
|++++++++++++++++++ | 72%
|
|+++++++++++++++++++ | 76%
|
|++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++ | 84%
|
|++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++ | 92%
|
|++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++| 100%
residuals2(fit_final, type="quantile") ### Gráfico de residuales y bandas de confianza
vdtest(fit_final) ### Supuesto de dispersión constante, se cumple##
## Score test for varying dispersion parameter
##
## Statistic = 0.13441
## degrees of freedom = 1
## p-value = 0.7139A partir de este resultado, diremos que DE se optimiza cuando se toman valores más pequeños de pH.
Modelo para el Porcentaje de Ácido Galacturónico (GA)
Generamos un modelo inicial para relacionar el porcentaje de GA con respecto a las variables explicativas del modelo.
mc3<- rsm(GA ~ SO(x1, x2, x3), data = MUCILAGO1)
summary(mc3)##
## Call:
## rsm(formula = GA ~ SO(x1, x2, x3), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.03774 1.92785 22.8429 4.611e-07 ***
## x1 -12.11611 1.36320 -8.8880 0.000113 ***
## x2 1.10665 1.36320 0.8118 0.447905
## x3 0.29589 1.36320 0.2171 0.835359
## x1:x2 0.53431 1.92785 0.2772 0.790963
## x1:x3 -2.00128 1.92785 -1.0381 0.339243
## x2:x3 -1.20131 1.92785 -0.6231 0.556127
## x1^2 0.36699 1.92785 0.1904 0.855304
## x2^2 1.61068 1.92785 0.8355 0.435449
## x3^2 2.28191 1.92785 1.1837 0.281314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9329, Adjusted R-squared: 0.8322
## F-statistic: 9.265 on 9 and 6 DF, p-value: 0.006775
##
## Analysis of Variance Table
##
## Response: GA
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 1184.90 394.97 26.5676 0.0007306
## TWI(x1, x2, x3) 3 22.93 7.64 0.5142 0.6873825
## PQ(x1, x2, x3) 3 31.74 10.58 0.7118 0.5796186
## Residuals 6 89.20 14.87
## Lack of fit 3 53.79 17.93 1.5188 0.3697992
## Pure error 3 35.41 11.80
##
## Stationary point of response surface:
## x1 x2 x3
## -83.479762 -0.191377 -36.721781
##
## Stationary point in original units:
## pH t_min inv_LS
## -39.239881 7.021558 -168.608903
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 3.01667788 1.30343436 -0.06053489
##
## $vectors
## [,1] [,2] [,3]
## x1 -0.3571018 -0.1604421 0.920182936
## x2 -0.4235173 0.9058653 -0.006411609
## x3 0.8325331 0.3920030 0.391436145Veamos la gráfica de la superficie de respuesta del modelo inicial:
par(mfrow=c(1,3))
contour (mc3, ~ x1 + x2 + x3, image=TRUE)
persp(mc3,~x1+x2+x3, col="lightblue",contours = list(z="top", col="orange"),theta = -145, phi = 35)
Eliminamos los efectos que generan ruido sobre el modelo y obtenemos el siguiente modelo
mc3.1<- rsm(GA ~ FO(x1), data = MUCILAGO1)
summary(mc3.1)##
## Call:
## rsm(formula = GA ~ FO(x1), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.16753 0.83017 55.612 < 2.2e-16 ***
## x1 -12.11611 1.17403 -10.320 6.314e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.8838, Adjusted R-squared: 0.8755
## F-statistic: 106.5 on 1 and 14 DF, p-value: 6.314e-08
##
## Analysis of Variance Table
##
## Response: GA
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1) 1 1174.40 1174.40 106.5039 6.314e-08
## Residuals 14 154.38 11.03
## Lack of fit 1 0.54 0.54 0.0455 0.8344
## Pure error 13 153.84 11.83
##
## Direction of steepest ascent (at radius 1):
## x1
## -1
##
## Corresponding increment in original units:
## pH
## -0.5shapiro.test(resid(mc3.1))##
## Shapiro-Wilk normality test
##
## data: resid(mc3.1)
## W = 0.86373, p-value = 0.02183bartlett.test(resid(mc3.1)~MUCILAGO1$x1)##
## Bartlett test of homogeneity of variances
##
## data: resid(mc3.1) by MUCILAGO1$x1
## Bartlett's K-squared = 0.022224, df = 2, p-value = 0.9889El modelo anterior no se distribuye normal, corregimos el modelo usando la transformación de la variable GA a escala log y obtenemos que el modelo se comporta bien y asimismo que cumple con todos los supuestos.
mc3.2<- rsm(log(GA) ~ FO(x1), data = MUCILAGO1)
summary(mc3.2)##
## Call:
## rsm(formula = log(GA) ~ FO(x1), data = MUCILAGO1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.811976 0.019791 192.6100 < 2.2e-16 ***
## x1 -0.269064 0.027989 -9.6132 1.519e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.8684, Adjusted R-squared: 0.859
## F-statistic: 92.41 on 1 and 14 DF, p-value: 1.519e-07
##
## Analysis of Variance Table
##
## Response: log(GA)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1) 1 0.57916 0.57916 92.4140 1.519e-07
## Residuals 14 0.08774 0.00627
## Lack of fit 1 0.00312 0.00312 0.4797 0.5007
## Pure error 13 0.08462 0.00651
##
## Direction of steepest ascent (at radius 1):
## x1
## -1
##
## Corresponding increment in original units:
## pH
## -0.5shapiro.test(resid(mc3.2))##
## Shapiro-Wilk normality test
##
## data: resid(mc3.2)
## W = 0.9226, p-value = 0.1857bartlett.test(resid(mc3.2)~MUCILAGO1$x1)##
## Bartlett test of homogeneity of variances
##
## data: resid(mc3.2) by MUCILAGO1$x1
## Bartlett's K-squared = 1.04, df = 2, p-value = 0.5945plot(MUCILAGO$pH,log(MUCILAGO1$GA))
Interpretación de resultados
Modelo Multivariado
Debemos tener presente que como el modelo mF es el escogido entonces la optimización conjunta de ambas variables (DE y GA) dependerá únicamente del pH.
MUCI<-MUCILAGO%>%group_by(pH)%>%summarise(Prom_DE=mean(DE), Prom_GA=mean(GA))
kable(x=MUCI,caption="Valores promedio para GA y DE para cada nivel de pH")| pH | Prom_DE | Prom_GA |
|---|---|---|
| 2.0 | 78.00000 | 58.46713 |
| 2.5 | 74.24375 | 45.98403 |
| 3.0 | 67.28250 | 34.23491 |
plot(MUCI$pH, MUCI$Prom_DE, ylim = c(35,80), xlab="Nivel de pH", ylab="%")
points(MUCI$pH,MUCI$Prom_GA, col="blue")
De acá podemos concluir que las variables DE y GA se optimizan
mutuamente cuando el pH es igual a 2. Esto hace que ambas tomen en
promedio valores más altos en dicho punto.
Adicionalmente las otras dos variables no tienen efecto sobre su resultado.
Modelos univariados
Para el rendimiento tenemos como significativos el pH y el tiempo de extracción, sin embargo, recomendamos tomar valores más pequeños para los niveles de pH y valores más altos del tiempo de extracción, debido a que no se puede decir con exactitud que sean los puntos donde se optimiza del todo el rendimiento. Por el momento, con los datos dados en el punto \(pH=2.0\) y \(t_{min}=10\) se optimiza el rendimiento.
Para el grado de esterificación y el ácido galacturónico sólo es significativo el pH, es importante tener cuidado porque el modelo del grado de esterificación no se comporta normal sino se ajusta con una distribución Gamma. En ambos casos las variables se optimizan en \(pH=2.0\). Se recomiendo evaluar valores más pequeños de pH o evaluar otras variables adicionales para observar un mejor comportamiento de ambas variables.
Referencias
- Montgomery, D. (2005). DISEÑO Y ANALISIS DE EXPERIMENTOS (2 Tra). Editorial Limusa.
- Myers, R. H., Montgomery, D. C. & Anderson-Cook, C. M. (2011). Response Surface Methodology: Process and Product Optimization Using Designed Experiments. Wiley.