Experimento de conversión [Diseño 2\(^4\) (2\(^k\))]
Kenneth Roy Cabrera Torres
Martes, 17 de abril de 2018
## Paquetes utilizados
library(readxl)
library(ggplot2)
library(pid)## Loading required package: png## Loading required package: FrF2## Loading required package: DoE.base## Loading required package: grid## Loading required package: conf.design##
## Attaching package: 'DoE.base'## The following objects are masked from 'package:stats':
##
## aov, lm## The following object is masked from 'package:graphics':
##
## plot.design## The following object is masked from 'package:base':
##
## lengthslibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Función auxiliar para hallar el error estándar de los resdiuales
# de un modelo "aov"
options(width=120)
errsdred<-function(modelo){
aovm<-anova(modelo)
return(sqrt(aovm$"Mean Sq"[length(aovm$"Mean Sq")]))
}En un proceso de desarrollo se quiere determinar cuáles de cuatro factores influyen más sobre la conversión de cierta materia prima en producto útil, medido en porcentaje.
Los factores a estudiar es la carga de catalizador (K lb), temperatura (Tp °C), presión (P psi) y concentración de un reactivo (C %).
Se aleatorizó el orden de ejecución de una sola corrida (1 réplica) y se halló el resultado de las 16 combianciones.
Datos tomados de libro Statistics for experimenters: design, innovation and discovery. Box G, Hunter J. y Hunter G. 2005. Tabla 5.10a.
Lectura y adecuación de la base de datos
library(readxl)
bd1 <- read_excel("conversion.xlsx")
bd1$K <- factor(bd1$K, labels = c("10 lb", "15 lb"))
bd1$Tp <- factor(bd1$Tp, labels = c("220 °C","240 °C"))
bd1$P <- factor(bd1$P, labels = c("50 psi", "80 psi"))
bd1$C <- factor(bd1$C, labels = c("10 %","12 %"))
bd1## # A tibble: 16 x 7
## OrdYates K Tp P C conversión OrdAleat
## <dbl> <fct> <fct> <fct> <fct> <dbl> <dbl>
## 1 1. 10 lb 220 °C 50 psi 10 % 70. 8.
## 2 2. 15 lb 220 °C 50 psi 10 % 60. 2.
## 3 3. 10 lb 240 °C 50 psi 10 % 89. 10.
## 4 4. 15 lb 240 °C 50 psi 10 % 81. 4.
## 5 5. 10 lb 220 °C 80 psi 10 % 69. 15.
## 6 6. 15 lb 220 °C 80 psi 10 % 62. 9.
## 7 7. 10 lb 240 °C 80 psi 10 % 88. 1.
## 8 8. 15 lb 240 °C 80 psi 10 % 81. 13.
## 9 9. 10 lb 220 °C 50 psi 12 % 60. 16.
## 10 10. 15 lb 220 °C 50 psi 12 % 49. 5.
## 11 11. 10 lb 240 °C 50 psi 12 % 88. 11.
## 12 12. 15 lb 240 °C 50 psi 12 % 82. 14.
## 13 13. 10 lb 220 °C 80 psi 12 % 60. 3.
## 14 14. 15 lb 220 °C 80 psi 12 % 52. 12.
## 15 15. 10 lb 240 °C 80 psi 12 % 86. 6.
## 16 16. 15 lb 240 °C 80 psi 12 % 79. 7.Análisis exploratorio
library(ggplot2)
g1 <- ggplot(bd1, aes(Tp, conversión, col = P))
g1 + geom_point() +
geom_line(aes(group = P)) +
facet_grid(C ~ K) +
scale_colour_discrete(name = "Presión") +
xlab("Temperatura") +
ylab("Conversión (%)")
library(ggplot2)
g2 <- ggplot(bd1, aes(P, conversión, col = Tp))
g2 + geom_point() +
geom_line(aes(group = Tp)) +
facet_grid(K ~ C) +
scale_colour_discrete(name = "Temperatura") +
xlab("Presión") +
ylab("Conversión (%)")
Modelo completo
El modelo propuesto completo sería \[ y_{ijklm} = \mu + \alpha_i + \beta_j + \gamma_k + \delta_l + \\ (\alpha \beta)_{ij} + (\alpha \gamma)_{ik} + (\alpha \delta)_{il} + (\beta \gamma)_{jk} + (\beta \delta)_{jl} + (\gamma \delta)_{kl} + \\ (\alpha \beta \gamma)_{ijk} + (\alpha \beta \delta)_{ijl} + (\beta \gamma \delta)_{jkl} + (\alpha \beta \delta)_{ijl} + \\ (\alpha \beta \gamma \delta)_{ijkl} + \epsilon_{ijklm} \] En este caso es un modelo saturado porque los grados de libertad, es igual al número de parámetros que se necesitan estimar.
Modelización
modelo1 <- aov(conversión ~ K + Tp + P + C +
K:Tp + K:P + K:C + Tp:P + Tp:C + P:C +
K:Tp:P + K:Tp:C + Tp:P:C + K:P:C +
K:Tp:P:C,
data = bd1)
summary(modelo1)## Df Sum Sq Mean Sq
## K 1 256.0 256.0
## Tp 1 2304.0 2304.0
## P 1 0.2 0.2
## C 1 121.0 121.0
## K:Tp 1 4.0 4.0
## K:P 1 2.2 2.2
## K:C 1 0.0 0.0
## Tp:P 1 6.3 6.3
## Tp:C 1 81.0 81.0
## P:C 1 0.3 0.3
## K:Tp:P 1 2.3 2.3
## K:Tp:C 1 1.0 1.0
## Tp:P:C 1 2.3 2.3
## K:P:C 1 0.2 0.2
## K:Tp:P:C 1 0.2 0.2Gráfica de Pareto
library(pid)
paretoPlot(modelo1)
efectos<-unlist(model.tables(modelo1))
efectos<-efectos[substring(names(efectos),1,6)=="tables"]
names(efectos)<-gsub("tables.","",names(efectos))
qqR<-qqnorm(efectos, xlim = c(-3, 3), ylim = c(-15,15), pch = 19, cex = 0.5)
qqline(efectos)
with(qqR,text(x,y,labels=names(efectos),cex=0.7,pos=3, srt=45))
modelo1 <- aov(conversión ~ K + Tp + P + C +
K:Tp + K:P + K:C + Tp:P + Tp:C + P:C +
K:Tp:P + K:Tp:C + Tp:P:C + K:P:C,
data = bd1)
summary(modelo1)## Df Sum Sq Mean Sq F value Pr(>F)
## K 1 256.0 256.0 1024 0.01989 *
## Tp 1 2304.0 2304.0 9216 0.00663 **
## P 1 0.2 0.2 1 0.50000
## C 1 121.0 121.0 484 0.02892 *
## K:Tp 1 4.0 4.0 16 0.15596
## K:P 1 2.2 2.2 9 0.20483
## K:C 1 0.0 0.0 0 1.00000
## Tp:P 1 6.3 6.3 25 0.12567
## Tp:C 1 81.0 81.0 324 0.03533 *
## P:C 1 0.3 0.3 1 0.50000
## K:Tp:P 1 2.3 2.3 9 0.20483
## K:Tp:C 1 1.0 1.0 4 0.29517
## Tp:P:C 1 2.3 2.3 9 0.20483
## K:P:C 1 0.2 0.2 1 0.50000
## Residuals 1 0.2 0.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo2 <- aov(conversión ~ K + Tp + P + C +
K:Tp +
Tp:P + Tp:C,
data = bd1)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## K 1 256.0 256.0 240.941 2.95e-07 ***
## Tp 1 2304.0 2304.0 2168.471 5.00e-11 ***
## P 1 0.2 0.2 0.235 0.6406
## C 1 121.0 121.0 113.882 5.21e-06 ***
## K:Tp 1 4.0 4.0 3.765 0.0883 .
## Tp:P 1 6.3 6.3 5.882 0.0415 *
## Tp:C 1 81.0 81.0 76.235 2.31e-05 ***
## Residuals 8 8.5 1.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo3 <- aov(conversión ~ K + Tp + C +
Tp:C,
data = bd1)
summary(modelo3)## Df Sum Sq Mean Sq F value Pr(>F)
## K 1 256 256.0 148.21 1.00e-07 ***
## Tp 1 2304 2304.0 1333.89 7.81e-13 ***
## C 1 121 121.0 70.05 4.24e-06 ***
## Tp:C 1 81 81.0 46.90 2.77e-05 ***
## Residuals 11 19 1.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1paretoPlot(modelo3)
Diagnóstico del modelo
Normalidad.
shapiro.test(rstandard(modelo3))##
## Shapiro-Wilk normality test
##
## data: rstandard(modelo3)
## W = 0.93303, p-value = 0.2722Homocedasticidad.
bartlett.test(conversión ~ K, data=bd1)##
## Bartlett test of homogeneity of variances
##
## data: conversión by K
## Bartlett's K-squared = 0.048882, df = 1, p-value = 0.825bartlett.test(conversión ~ Tp, data=bd1)##
## Bartlett test of homogeneity of variances
##
## data: conversión by Tp
## Bartlett's K-squared = 2.3461, df = 1, p-value = 0.1256bartlett.test(conversión ~ C, data=bd1)##
## Bartlett test of homogeneity of variances
##
## data: conversión by C
## Bartlett's K-squared = 0.75001, df = 1, p-value = 0.3865bartlett.test(conversión ~ paste(Tp,C), data=bd1)##
## Bartlett test of homogeneity of variances
##
## data: conversión by paste(Tp, C)
## Bartlett's K-squared = 0.34548, df = 3, p-value = 0.9513Independencia.
library(lmtest)
dwtest(conversión ~ Tp + K + C + Tp:C, order.by = ~ OrdAleat, data = bd1,
alternative = "two.sided")##
## Durbin-Watson test
##
## data: conversión ~ Tp + K + C + Tp:C
## DW = 1.3816, p-value = 0.2694
## alternative hypothesis: true autocorrelation is not 0Inferencia
model.tables(modelo3, se = TRUE)## Tables of effects
##
## K
## K
## 10 lb 15 lb
## 4 -4
##
## Tp
## Tp
## 220 °C 240 °C
## -12 12
##
## C
## C
## 10 % 12 %
## 2.75 -2.75
##
## Tp:C
## C
## Tp 10 % 12 %
## 220 °C 2.25 -2.25
## 240 °C -2.25 2.25
##
## Standard errors of effects
## K Tp C Tp:C
## 0.4647 0.4647 0.4647 0.6571
## replic. 8 8 8 4model.tables(modelo3, type = "means", se = TRUE)## Tables of means
## Grand mean
##
## 72.25
##
## K
## K
## 10 lb 15 lb
## 76.25 68.25
##
## Tp
## Tp
## 220 °C 240 °C
## 60.25 84.25
##
## C
## C
## 10 % 12 %
## 75.0 69.5
##
## Tp:C
## C
## Tp 10 % 12 %
## 220 °C 65.25 55.25
## 240 °C 84.75 83.75
##
## Standard errors for differences of means
## K Tp C Tp:C
## 0.6571 0.6571 0.6571 0.9293
## replic. 8 8 8 4Gráficas del modelo final
bd1 <- cbind(bd1, predict(modelo3, interval = "confidence"))
library(ggplot2)
g1 <- ggplot(bd1, aes(Tp,conversión, col = C))
g1 + geom_point() +
facet_grid(. ~ K) +
scale_colour_discrete(name = "Concentración") +
geom_errorbar(aes(ymax = upr, ymin = lwr, col = C), width = 0.25) +
stat_summary(aes(Tp, fit, group = C, col = C), fun.y = mean,
geom = "line", size = 1) +
ylab("Conversión (Porcentaje)") +
xlab("Temperatura")
Estimación de la desviación estándar de error aleatorio
errsdred(modelo3)## [1] 1.314257