Taller de Práctica: Diseños Factoriales 3k

Maestría en Investigación Operativa y Estadística

Ejercicio 2: Mejorando la Resistencia de un Material Compuesto

Contexto

Una empresa de materiales desea mejorar la resistencia a la tracción de un nuevo material compuesto. Tres factores se consideran críticos:

Factor A: Proporción de fibra de vidrio (%)

Niveles: - Bajo (-1): 10% - Medio (0): 20% - Alto (+1): 30%

Factor B: Tipo de resina

Niveles: - Bajo (-1): Resina A - Medio (0): Resina B - Alto (+1): Resina C

Factor C: Temperatura de curado (°C)

Niveles: - Bajo (-1): 100°C - Medio (0): 125°C - Alto (+1): 150°C

Se realiza un diseño factorial completo (3^k) con 2 réplicas para cada tratamiento, obteniendo 54 observaciones en total.

Lectura de datos

datos_2 <- read_xlsx("datos_taller3.xlsx", sheet = "ejercicio2") %>% clean_names()
#ftable(datos)

Reetiquetado de variables

# A: fibra_de_vidrio
# B: b_resina
# C: c_temperatura

colnames(datos_2) <- c('tratamiento', 'replica', 'A', 'B', 'C', 'resistencia_m_pa')

datos_2 <- datos_2 %>% mutate(a_fibra_de_vidrio = case_when(A == -1 ~ 0.1,
                                            A == 0 ~ 0.2,
                                            TRUE ~ 0.3),
                              b_resina = case_when(B == -1 ~ "Resina_A",
                                            B == 0 ~ "Resina_B",
                                            TRUE ~ "Resina_C"),
                              c_temperatura = case_when(C == -1 ~ 100,
                                            C == 0 ~ 125,
                                            TRUE ~ 150))

datos_2

## # A tibble: 54 × 9
##    tratamiento replica     A     B     C resistencia_m_pa a_fibra_de_vidrio
##          <dbl>   <dbl> <dbl> <dbl> <dbl>            <dbl>             <dbl>
##  1           1       1    -1    -1    -1             49.4               0.1
##  2           1       2    -1    -1    -1             51.2               0.1
##  3           2       1    -1    -1     0             55.7               0.1
##  4           2       2    -1    -1     0             56.9               0.1
##  5           3       1    -1    -1     1             54.7               0.1
##  6           3       2    -1    -1     1             54.3               0.1
##  7           4       1    -1     0    -1             58.6               0.1
##  8           4       2    -1     0    -1             58.9               0.1
##  9           5       1    -1     0     0             62.7               0.1
## 10           5       2    -1     0     0             62.8               0.1
## # ℹ 44 more rows
## # ℹ 2 more variables: b_resina <chr>, c_temperatura <dbl>

Modelo lineal

modelo <- lm(resistencia_m_pa ~ A + B + C +I(A^2) + I(C^2)+ A:B +A:C + B:C, data=datos_2)
summary(modelo)

## 
## Call:
## lm(formula = resistencia_m_pa ~ A + B + C + I(A^2) + I(C^2) + 
##     A:B + A:C + B:C, data = datos_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.81806 -0.64306  0.07639  0.78403  2.01806 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 80.47778    0.33470 240.447  < 2e-16 ***
## A           15.30000    0.18332  83.459  < 2e-16 ***
## B            9.90833    0.18332  54.048  < 2e-16 ***
## C            1.93611    0.18332  10.561 9.14e-14 ***
## I(A^2)      -0.41667    0.31752  -1.312    0.196    
## I(C^2)      -3.09167    0.31752  -9.737 1.19e-12 ***
## A:B          0.02917    0.22452   0.130    0.897    
## A:C          0.18750    0.22452   0.835    0.408    
## B:C         -0.18333    0.22452  -0.817    0.418    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.1 on 45 degrees of freedom
## Multiple R-squared:  0.9956, Adjusted R-squared:  0.9948 
## F-statistic:  1262 on 8 and 45 DF,  p-value: < 2.2e-16

Ajuste de la ANOVA

anova(modelo)

## Analysis of Variance Table
## 
## Response: resistencia_m_pa
##           Df Sum Sq Mean Sq   F value    Pr(>F)    
## A          1 8427.2  8427.2 6965.4429 < 2.2e-16 ***
## B          1 3534.3  3534.3 2921.2390 < 2.2e-16 ***
## C          1  134.9   134.9  111.5389 9.144e-14 ***
## I(A^2)     1    2.1     2.1    1.7220    0.1961    
## I(C^2)     1  114.7   114.7   94.8047 1.191e-12 ***
## A:B        1    0.0     0.0    0.0169    0.8972    
## A:C        1    0.8     0.8    0.6974    0.4081    
## B:C        1    0.8     0.8    0.6667    0.4185    
## Residuals 45   54.4     1.2                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interacción A:B

library(ggplot2)

ggplot(datos_2, aes(x = factor(A), y = resistencia_m_pa, color = factor(B))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = B)) +
  labs(title = "Interacción entre Temperatura de curado (°C) y resistencia_m_pa",
       x = " Temperatura de curado (°C)(codificada)",
       y = "resistencia_m_pa",
       color = "Temperatura de curado (°C)") +
  theme_minimal()

Interacción A:C

library(ggplot2)

ggplot(datos_2, aes(x = factor(A), y = resistencia_m_pa, color = factor(C))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = C)) +
  labs(title = "Interacción entre Proporción de fibra de vidrio  y resistencia_m_pa",
       x = " Proporción de fibra de vidrio (codificada)",
       y = "resistencia_m_pa",
       color = "Proporción de fibra de vidrio") +
  theme_minimal()

## Interacción B:C

library(ggplot2)

ggplot(datos_2, aes(x = factor(B), y = resistencia_m_pa, color = factor(C))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = C)) +
  labs(title = "Interacción entre tipo de resina  y resistencia_m_pa",
       x = " Tipo de resina (codificada)",
       y = "resistencia_m_pa",
       color = "Proporción de fibra de vidrio") +
  theme_minimal()

## Superficie de respuesta

library(rsm)

# Ajustar el modelo para rsm
modelo_rsm <- rsm(resistencia_m_pa ~ SO(A, B, C), data = datos_2)

# Superficie de respuesta para A y B, con C en nivel central (0)
par(mfrow = c(1, 2))

# Contorno
contour(modelo_rsm, ~ A + C, at = list(B = 0), image = TRUE,
        main = "Contorno: resistencia_m_pa vs A y C (B=0)")

# Perspectiva 3D
persp(modelo_rsm, ~ A + C, at = list(B = 0), col = "lightblue",
      main = "Superficie: resistencia_m_pa vs A y C (B=0)")

Encontrar el punto óptimo

optimo_2 <- optim(
  par = c(A = 0, B = 0, C = 0), # Valores iniciales
  fn = function(x) {
    predict(modelo, newdata = data.frame(
      A = x[1],
      B = x[2],
      C = x[3]
    ))
  },
  control = list(fnscale = -1) # Maximizar la función
)

# Valores óptimos codificados
optimo_2$par

##             A             B             C 
##  1.579577e+29  1.617046e+30 -1.918718e+28

Conclusiones Caso 2

No es posible encontrar un valor óptimo dada la naturaleza (forma) de la superficie de respuesta. Existe sólo un efecto cuadrático significativo que le da cierta curvatura al “plano”.
No existe evidencia estadística con la cual poder establecer interacción entre la proporción de fibra de vidrio, el tipo de resina y la temperatura. En ese sentido, puede decirse que son independientes.
La proporción de fibra de vidrio tiene una mayor efecto sobre la resistencia de la resina que el tipo de resina (54.42% mayor), y el tipo de resina a su vez, tiene un mayor efecto sobre la resistencia de la resina que la temperatura (411.78% mayor).

Ejercicio 3: Desaroollo de un proceso químico

Contexto

Un ingeniero químico está optimizando un proceso de síntesis y necesita estudiar el efecto de tres factores en el rendimiento del producto:

Factor A: pH de la solución

Niveles: - Bajo (-1): 5 - Medio (0): 7 - Alto (+1): 9

Factor B: Concentración del reactivo (%)

Niveles: - Bajo (-1): 10% - Medio (0): 15% - Alto (+1): 20%

Factor C: Tiempo de reacción (min)

Niveles: - Bajo (-1): 30 min - Medio (0): 45 min - Alto (+1): 60 min

Se sospecha que existen interacciones significativas, especialmente entre el pH y la concentración del reactivo.

Lectura de datos

datos_3 <- read_xlsx("datos_taller3.xlsx", sheet = "ejercicio3") %>% clean_names()
#ftable(datos)

Reetiquetado de las variables

# A: a_p_h
# B: b_concentracion
# C: c_tiempo

colnames(datos_3) <- c('tratamiento', 'A', 'B', 'C', 'rendimiento_percent')

Asignación de valores iniciales

datos_3 <- datos_3 %>% mutate(a_p_h = case_when(A == -1 ~ 5,
                                            A == 0 ~ 7,
                                            TRUE ~ 9),
                              b_concentracion = case_when(B == -1 ~ 0.1,
                                            B == 0 ~ 0.15,
                                            TRUE ~ 0.2),
                              c_tiempo = case_when(C == -1 ~ 30,
                                            C == 0 ~ 45,
                                            TRUE ~ 60))

datos_3

## # A tibble: 27 × 8
##    tratamiento     A     B     C rendimiento_percent a_p_h b_concentracion
##          <dbl> <dbl> <dbl> <dbl>               <dbl> <dbl>           <dbl>
##  1           1    -1    -1    -1                60.6     5            0.1 
##  2           2    -1    -1     0                62.3     5            0.1 
##  3           3    -1    -1     1                63.4     5            0.1 
##  4           4    -1     0    -1                68.5     5            0.15
##  5           5    -1     0     0                69.3     5            0.15
##  6           6    -1     0     1                72.3     5            0.15
##  7           7    -1     1    -1                80.9     5            0.2 
##  8           8    -1     1     0                82.1     5            0.2 
##  9           9    -1     1     1                84.6     5            0.2 
## 10          10     0    -1    -1                65.6     7            0.1 
## # ℹ 17 more rows
## # ℹ 1 more variable: c_tiempo <dbl>

Modelo Lineal

modelo_3 <- lm(rendimiento_percent ~ A + B + C +I(A^2) + I(C^2) + I(C^2)+ A:B +A:C + B:C, data=datos_3)
summary(modelo)

## 
## Call:
## lm(formula = resistencia_m_pa ~ A + B + C + I(A^2) + I(C^2) + 
##     A:B + A:C + B:C, data = datos_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.81806 -0.64306  0.07639  0.78403  2.01806 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 80.47778    0.33470 240.447  < 2e-16 ***
## A           15.30000    0.18332  83.459  < 2e-16 ***
## B            9.90833    0.18332  54.048  < 2e-16 ***
## C            1.93611    0.18332  10.561 9.14e-14 ***
## I(A^2)      -0.41667    0.31752  -1.312    0.196    
## I(C^2)      -3.09167    0.31752  -9.737 1.19e-12 ***
## A:B          0.02917    0.22452   0.130    0.897    
## A:C          0.18750    0.22452   0.835    0.408    
## B:C         -0.18333    0.22452  -0.817    0.418    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.1 on 45 degrees of freedom
## Multiple R-squared:  0.9956, Adjusted R-squared:  0.9948 
## F-statistic:  1262 on 8 and 45 DF,  p-value: < 2.2e-16

Ajuste de la ANOVA

anova(modelo_3)

## Analysis of Variance Table
## 
## Response: rendimiento_percent
##           Df  Sum Sq Mean Sq  F value    Pr(>F)    
## A          1  305.87  305.87  71.6269 1.089e-07 ***
## B          1 1208.68 1208.68 283.0429 1.862e-12 ***
## C          1   67.28   67.28  15.7553 0.0008993 ***
## I(A^2)     1    0.81    0.81   0.1889 0.6689984    
## I(C^2)     1    0.17    0.17   0.0390 0.8456047    
## A:B        1   15.87   15.87   3.7164 0.0698114 .  
## A:C        1    0.19    0.19   0.0439 0.8363785    
## B:C        1    0.04    0.04   0.0096 0.9231827    
## Residuals 18   76.87    4.27                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interacción A:B

library(ggplot2)

ggplot(datos_3, aes(x = factor(A), y = rendimiento_percent, color = factor(B))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = B)) +
  labs(title = "Interacción entre a_p_h y rendimiento_percent",
       x = " a_p_h (codificada)",
       y = "rendimiento_percent",
       color = "b_concentracion") +
  theme_minimal()

Interacción A:C

library(ggplot2)

ggplot(datos_3, aes(x = factor(A), y = rendimiento_percent, color = factor(C))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = C)) +
  labs(title = "Interacción entre a_p_h y rendimiento_percent",
       x = " a_p_h (codificada)",
       y = "rendimiento_percent",
       color = "tiempo de curado") +
  theme_minimal()

Interacción B:C

library(ggplot2)

ggplot(datos_3, aes(x = factor(B), y = rendimiento_percent, color = factor(C))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = C)) +
  labs(title = "Interacción entre b_concentracion y rendimiento_percent",
       x = " b_concentracion (codificada)",
       y = "rendimiento_percent",
       color = "tiempo de curado") +
  theme_minimal()

Superficie de respuesta

library(rsm)

# Ajustar el modelo para rsm
modelo_rsm <- rsm(rendimiento_percent ~ SO(A, B, C), data = datos_3)

# Superficie de respuesta para A y B, con C en nivel central (0)
par(mfrow = c(1, 2))

# Contorno
contour(modelo_rsm, ~ A + B + C, at = list(A = 1), image = TRUE,
        main = "Contorno: rendimiento_percent vs A y C (B=0)")

# Perspectiva 3D
persp(modelo_rsm, ~ A + B + C, at = list(A = 1), col = "lightblue",
      main = "Superficie: resistencia_m_pa vs A y C (B=0)")

Encuentro del punto óptimo

optimo_3 <- optim(
  par = c(A = 0, B = 0, C = 0), # Valores iniciales
  fn = function(x) {
    predict(modelo_3, newdata = data.frame(
      A = x[1],
      B = x[2],
      C = x[3]
    ))
  },
  control = list(fnscale = -1) # Maximizar la función
)

# Valores óptimos codificados
optimo_3$par

##             A             B             C 
## -1.777758e+41  3.813311e+41  9.907990e+40

Conclusiones caso 3

No es posible encontrar un valor óptimo dada la naturaleza (forma) de la superficie de respuesta. Los efectos cuadrados no son significativos.
No existe evidencia estadística con la cual poder establecer interacción entre el pH de la solución, la concentración del reactivo y el tiempo de reacción. En ese sentido, puede decirse que son independientes.
La concentración del reactivo tiene una mayor efecto sobre el rendimiento del producto que el pH de la solución (98.79% mayor), y el pH de la solución a su vez, tiene un mayor efecto sobre el rendimiento del producto que el tiempo de reacción (113.24% mayor).