Diseño y Análisis de Experimentos

Resultados Trabajo #2

Autor/a

Afiliación

Jason E. Barrios Mosquera

Acesco Colombia SAS, Malambo

Homework 2

1. Reproducir los resultados de Example 6.1 y Example 6.2.
2. Resolver Problem 6.1, Problem 6.5 y Problem 6.15.

de la \(9^a\) edición de Design of Experiments por Douglas Montgomery. Para más información, ver aquí.

Fecha de entrega: 30 de Agosto.

Ejemplo 6.1

Lectura de datos

Los datos con los que trabajaremos son los siguientes:

## read the data set
datos <- read.table('C:/Users/FBA/OneDrive - Acesco/0. Personales/2. Profesional/20240416 - Maestria en Analítica de Datos/20240726 Diseño de Experimentos/Homework 2/Datos Example 6.1.txt', header = TRUE) 

datos$A<-as.factor(datos$A)
datos$B<-as.factor(datos$B)
datos$C<-as.factor(datos$C)
datos$Run<-as.character(datos$Run)

## see first 5 rows of the data
head(datos, 20)

   Run A B C Yield
1    1 - - -   550
2    2 + - -   669
3    3 - + -   633
4    4 + + -   642
5    5 - - +  1037
6    6 + - +   749
7    7 - + +  1075
8    8 + + +   729
9    9 - - -   604
10  10 + - -   650
11  11 - + -   601
12  12 + + -   635
13  13 - - +  1052
14  14 + - +   868
15  15 - + +  1063
16  16 + + +   860

El número total de experimentos es:

## número de experimentos
NROW(datos)

[1] 16

Para obtener el número de réplicas para cada nivel de A, B y C hacemos:

## número de réplicas
with(datos, tapply(Yield, list(A, B, C), length))

, , -

  - +
- 2 2
+ 2 2

, , +

  - +
- 2 2
+ 2 2

Análisis gráfico

Ahora, procedemos a realizar un análisis gráfico de los resultados del experimento.

Inicialmente, los boxplots para los efectos principales son:

# boxplots for main effects
boxplot(Yield ~ A*B*C, datos, las = 1, col = 2:4)

El gráfico de interacción puede obtenerse haciendo:

# interaction plot
with(datos, interaction.plot(A, B, Yield, las = 1))

with(datos, interaction.plot(B, C, Yield, las = 1))

with(datos, interaction.plot(A, C, Yield, las = 1))

Graficamente parece que los factores A y C si afectan significativamente los resultados de la variable respuesta y parece existir interacción significativa. Adicionalemente parece existir interacción significativa entre A y B. Esto se validará a través de la tabla ANOVA.

Construcción de la tabla ANOVA

Inicialmente consideraremos la tabla ANOVA con interacción

## full ANOVA model
fit1 <- aov(Yield ~ A * B * C, datos)
summary(fit1)

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1  41311   41311  18.339 0.002679 ** 
B            1    218     218   0.097 0.763911    
C            1 374850  374850 166.411 1.23e-06 ***
A:B          1   2475    2475   1.099 0.325168    
A:C          1  94403   94403  41.909 0.000193 ***
B:C          1     18      18   0.008 0.930849    
A:B:C        1    127     127   0.056 0.818586    
Residuals    8  18020    2253                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso, no tiene sentido considerar la interacción entre A:B ni, B:C, ni A:B:C, ni B como un factor adicional que explique la variable Yield.

Al removemos los términos no significativos para \(\alpha=0.05\)

## reduced ANOVA model
fit <- aov(Yield ~ A +C+A:C, datos)
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1  41311   41311   23.77 0.000382 ***
C            1 374850  374850  215.66 4.95e-09 ***
A:C          1  94403   94403   54.31 8.62e-06 ***
Residuals   12  20858    1738                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Por lo tanto, concluimos que ambos A y C y A:C factores tienen un efecto sobre la variable Yield: Etch rate.

Análisis del experimento utilizando RLM

Una vez detectamos que los factores de interés tienen un efecto sobre la variable respuesta Yield, es importante cuantificar dicho efecto.

Con el fin de utilizar RLM debemos recodificar los factores A y B como

\[ \begin{eqnarray} x_1 = \frac{A-(A^- + A^+)/2}{(A^+ - A^-)/2} \nonumber \\ x_2 = \frac{C-(C^- + C^+)/2}{(C^+ - C^-)/2} \nonumber \end{eqnarray} \]

Para recodificar las variables en R hacemos:

## primero codificamos los factores para que estén en el rango -1 y +1 
## (ver # arriba)
coded <- function(x) ifelse(x == x[1], -1, 1)
datos <- within(datos, {
    coded_A = coded(A)
    coded_C = coded(C) 
    coded_B = coded(B)
    })
head(datos, 10)

   Run A B C Yield coded_B coded_C coded_A
1    1 - - -   550      -1      -1      -1
2    2 + - -   669      -1      -1       1
3    3 - + -   633       1      -1      -1
4    4 + + -   642       1      -1       1
5    5 - - +  1037      -1       1      -1
6    6 + - +   749      -1       1       1
7    7 - + +  1075       1       1      -1
8    8 + + +   729       1       1       1
9    9 - - -   604      -1      -1      -1
10  10 + - -   650      -1      -1       1

Posteriormente estimamos el modelo de RLM haciendo:

# now fit the regression excluding the interaction term
fitlm1 <- lm(Yield ~ coded(A)*coded(B)*coded(C), datos)
fitlm <- lm(Yield ~ coded(A) + coded(C) + coded(A):coded(C), datos)
summary(fitlm)

Call:
lm(formula = Yield ~ coded(A) + coded(C) + coded(A):coded(C), 
    data = datos)

Residuals:
   Min     1Q Median     3Q    Max 
-72.50 -15.44   2.50  18.69  66.50 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         776.06      10.42  74.458  < 2e-16 ***
coded(A)            -50.81      10.42  -4.875 0.000382 ***
coded(C)            153.06      10.42  14.685 4.95e-09 ***
coded(A):coded(C)   -76.81      10.42  -7.370 8.62e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 41.69 on 12 degrees of freedom
Multiple R-squared:  0.9608,    Adjusted R-squared:  0.9509 
F-statistic: 97.91 on 3 and 12 DF,  p-value: 1.054e-08

Observe que los para ́metros del modelo ajustado son

\(\hat{\mathbf{\beta}} = (\) 776.0625, -50.8125, 153.0625, -76.8125\()\)

que corresponden a las cantidades anteriormente calculadas a partir de los efectos principales.

Intervalos de confianza

confint(fitlm1)

                                2.5 %    97.5 %
(Intercept)                 748.70109 803.42391
coded(A)                    -78.17391 -23.45109
coded(B)                    -23.67391  31.04891
coded(C)                    125.70109 180.42391
coded(A):coded(B)           -39.79891  14.92391
coded(A):coded(C)          -104.17391 -49.45109
coded(B):coded(C)           -28.42391  26.29891
coded(A):coded(B):coded(C)  -24.54891  30.17391

Los intervalos de confianza ratifican la desición anterior. Los únicos factores cuyo intervalo no incluy el cero son A y C y su interacción A:C

Validación de supuestos

En el DOE \(2^k\) deben validarse los mismos supuestos sobre el error. En R dicha validación se realiza sobre los residuales del modelo ajustado, que en nuestro corresponde al objeto fitlm.

Para ello utilizaremos la función validar_supuestos() sobre el objeto fitlm. Primero cargamos la función utilizando source()

## cargar la función
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')

Esta función fue creada por el profesor Jorge Vélez, PhD 
<https://jorgeivanvelez.netlify.app/> para la validación 
de los supuestos básicos en RLM. 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a <jvelezv@uninorte.edu.co> 
 
Última modificación:  Agosto 13, 2024 
 

Finalmente aplicamos la función:

## validación de supuestos
validar_supuestos(fitlm)

      normalidad homocedasticidad    independencia 
        "cumple"         "cumple"         "cumple" 

Conclusión: Los resultados indican que los errores del modelo ajustado son independientes, siguen una distribución Normal y tienen varianza constante.

La conclusión global de este análisis es que los factores A, C y su interacción tienen un efecto importante sobre la variable :Etch rate y . Como los residuales del modelo ajustado cumplen con todos los supuestos, las conclusiones que se deriven de él son válidas para todo el proceso.

Superficie de Respuesta

A partir del modelo de RLM ajustado, contenido en el objeto fitlm, es posible utilizar grid-search para construir la superficie de respuesta asociada y determinar las condiciones de operación óptimas que garantizan un valor máximo de Yield: Etch rate.

Inicialmente creamos la rejilla de evaluación:

# grid-based optimization
ngrid <- 30
yieldLm <- lm(Yield ~ coded_A + coded_C + coded_A:coded_C, data = datos)
Agrid <- Cgrid <- seq(from = -1, to = 1, length = ngrid)
yhat <- predict(yieldLm, expand.grid(coded_A = Agrid, coded_C = Cgrid))
yhat <- matrix(yhat, length(Agrid), length(Cgrid))

y posteriormente procedemos a graficar la superficie y los contornos:

## disponibilidad del paquete plot3D 
if (!require(plot3D)) install.packages("plot3D")
require(plot3D)

## rearranging the window
par(mfrow = c(1, 2))

## 3D plot (same as Fig 6.3a)
persp3D(((1.2+0.8)/2)+(1.2-0.8)*Agrid,((325+275)/2)+(325-275)*Cgrid, yhat, 
        ticktype = "detailed", 
        colkey = FALSE, phi = 30, theta = -40, bty = "b2", 
        space = 0.15, d = 1.5, 
        xlab = "\n A: Gap", 
        ylab = "\n C: Power", 
        zlab = "\n\n Etch rate", contour = TRUE, 
        zlim = c(500, 1100))

## 2D plot (same as Fig 6.3b)
contour2D(yhat, ((1.2+0.8)/2)+(1.2-0.8)*Agrid,((325+275)/2)+(325-275)*Cgrid, 
        xlab = "\n A: Gap", 
        ylab = "\n C: Power",  colkey = FALSE, las = 1)

# Definir los rangos de A y C con 100 puntos entre -1 y 1
A_values <- seq(-1, 1, length.out = 100)
C_values <- seq(-1, 1, length.out = 100)

# Crear una cuadrícula de valores para A y C
grid <- expand.grid(A = A_values, C = C_values)

# Calcular el término de interacción A:C
grid$AC <- grid$A * grid$C

# Calcular Yield para cada combinación de A y C
grid$Yield <- summary(fitlm)[[4]][1] + summary(fitlm)[[4]][2] * grid$A + summary(fitlm)[[4]][3] * grid$C + summary(fitlm)[[4]][4] * grid$AC

# Convertir los valores de Yield en una matriz para la visualización
Yield_matrix <- matrix(grid$Yield, nrow = 100, ncol = 100)

# Crear la visualización de superficie
library(plotly)

Cargando paquete requerido: ggplot2

Adjuntando el paquete: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

fig <- plot_ly(z = ~Yield_matrix, x = ~A_values, y = ~C_values)
fig <- fig %>% add_surface()
fig <- fig %>% layout(scene = list(
  xaxis = list(title = "A: Gap"),
  yaxis = list(title = "C: Power"),
  zaxis = list(title = "Etch rate")
))
# Mostrar la gráfica
fig

## availability of knitr
if(!require(plotly)) install.packages('plotly')
require(plotly)

fig <- plot_ly(z = ~Yield_matrix)
fig <- fig %>% add_surface()

fig

Ejemplo 6.2

Lectura de datos

Los datos2 con los que trabajaremos son los siguientes:

## read the data set
datos2 <- read.table('C:/Users/FBA/OneDrive - Acesco/0. Personales/2. Profesional/20240416 - Maestria en Analítica de Datos/20240726 Diseño de Experimentos/Homework 2/Datos Example 6.2.txt', header = TRUE) 

datos2$A<-as.factor(datos2$A)
datos2$B<-as.factor(datos2$B)
datos2$C<-as.factor(datos2$C)
datos2$D<-as.factor(datos2$D)
datos2$RunLabel<-as.character(datos2$RunLabel)

## see first 5 rows of the data
head(datos2, 20)

   RunNumber A B C D RunLabel FiltrationRate
1          1 − − − −      (1)             45
2          2 + − − −        a             71
3          3 − + − −        b             48
4          4 + + − −       ab             65
5          5 − − + −        c             68
6          6 + − + −       ac             60
7          7 − + + −       bc             80
8          8 + + + −      abc             65
9          9 − − − +        d             43
10        10 + − − +       ad            100
11        11 − + − +       bd             45
12        12 + + − +      abd            104
13        13 − − + +       cd             75
14        14 + − + +      acd             86
15        15 − + + +      bcd             70
16        16 + + + +     abcd             96

El número total de experimentos es:

## número de experimentos
NROW(datos2)

[1] 16

Para obtener el número de réplicas para cada nivel de A, B y C hacemos:

## número de réplicas
with(datos2, tapply(FiltrationRate, list(A, B, C,D), length))

, , −, −

  − +
− 1 1
+ 1 1

, , +, −

  − +
− 1 1
+ 1 1

, , −, +

  − +
− 1 1
+ 1 1

, , +, +

  − +
− 1 1
+ 1 1

Análisis gráfico

Ahora, procedemos a realizar un análisis gráfico de los resultados del experimento.

Inicialmente, los boxplots para los efectos principales son:

# boxplots for main effects
boxplot(FiltrationRate ~ A*B*C*D, datos2, las = 1, col = 2:4)

El gráfico de interacción por pares puede obtenerse haciendo:

# interaction plot
with(datos2, interaction.plot(A, B, FiltrationRate, las = 1))

with(datos2, interaction.plot(B, C, FiltrationRate, las = 1))

with(datos2, interaction.plot(A, C, FiltrationRate, las = 1))

with(datos2, interaction.plot(D, A, FiltrationRate, las = 1))

with(datos2, interaction.plot(D, B, FiltrationRate, las = 1))

with(datos2, interaction.plot(D, C, FiltrationRate, las = 1))

with(datos2, plot(A,FiltrationRate, las = 1))

with(datos2, plot(B, FiltrationRate, las = 1))

with(datos2, plot(C, FiltrationRate, las = 1))

with(datos2, plot(D,  FiltrationRate, las = 1))

Construcción de la tabla ANOVA

Inicialmente consideraremos la tabla ANOVA con interacción

# Ajustar el modelo ANOVA
fit <- aov(FiltrationRate ~ A * B * C * D, data = datos2)

# Obtener el resumen del modelo
f <- summary(fit)

# Extraer la tabla ANOVA del resumen
aov_table <- f[[1]]

# Calcular la suma total de cuadrados
total_sum_sq <- sum(aov_table$`Sum Sq`)

# Calcular la contribución porcentual
aov_table$PercentContribution <- aov_table$`Sum Sq` / total_sum_sq * 100
aov_table$"***" <- ifelse(aov_table$PercentContribution > 1, 1, 0)

# Imprimir la tabla modificada
print(aov_table)

            Df  Sum Sq Mean Sq PercentContribution ***
A            1 1870.56 1870.56              32.640   1
B            1   39.06   39.06               0.682   0
C            1  390.06  390.06               6.806   1
D            1  855.56  855.56              14.929   1
A:B          1    0.06    0.06               0.001   0
A:C          1 1314.06 1314.06              22.929   1
B:C          1   22.56   22.56               0.394   0
A:D          1 1105.56 1105.56              19.291   1
B:D          1    0.56    0.56               0.010   0
C:D          1    5.06    5.06               0.088   0
A:B:C        1   14.06   14.06               0.245   0
A:B:D        1   68.06   68.06               1.188   1
A:C:D        1   10.56   10.56               0.184   0
B:C:D        1   27.56   27.56               0.481   0
A:B:C:D      1    7.56    7.56               0.132   0

*** = 1Factores que son significantes respecto a la proporción del error que lograr explicar. En este caso, no tiene sentido considerar B, A:B,B:C, B:D,C:D, A:B:C, A:C:D, B:C:D ni A:B:C:D como factores que que expliquen la respuesta FiltrationRate.

Al removemos los términos no significativos. Se incluye B por convneción de jerarquía,dada que su interacción con A y D parece ser significativa

## reduced ANOVA model
fit <- aov(FiltrationRate ~ A+B+C+D+A:C+A:D+A:B:D, datos2)
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1 1870.6  1870.6 128.451 2.83e-05 ***
B            1   39.1    39.1   2.682 0.152576    
C            1  390.1   390.1  26.785 0.002063 ** 
D            1  855.6   855.6  58.751 0.000258 ***
A:C          1 1314.1  1314.1  90.236 7.76e-05 ***
A:D          1 1105.6  1105.6  75.918 0.000126 ***
A:B:D        3   68.7    22.9   1.572 0.291182    
Residuals    6   87.4    14.6                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Por lo tanto, concluimos que ambos By A:B:D factores no tienen un efecto sobre la variable FiltrationRate, para \(\alpha=0.05\)

## reduced ANOVA model
fit <- aov(FiltrationRate ~ A+C+D+A:C+A:D, datos2)
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1 1870.6  1870.6   95.86 1.93e-06 ***
C            1  390.1   390.1   19.99   0.0012 ** 
D            1  855.6   855.6   43.85 5.92e-05 ***
A:C          1 1314.1  1314.1   67.34 9.41e-06 ***
A:D          1 1105.6  1105.6   56.66 2.00e-05 ***
Residuals   10  195.1    19.5                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis del experimento utilizando RLM

Una vez detectamos que los factores de interés tienen un efecto sobre la variable respuesta FiltrationRate, es importante cuantificar dicho efecto.

Con el fin de utilizar RLM debemos recodificar los factores A y B como

\[ \begin{eqnarray} x_1 = \frac{A-(A^- + A^+)/2}{(A^+ - A^-)/2} \nonumber \\ x_2 = \frac{C-(C^- + C^+)/2}{(C^+ - C^-)/2} \nonumber \\ x_3 = \frac{D-(D^- + D^+)/2}{(D^+ - D^-)/2} \nonumber\\ \end{eqnarray} \] Para recodificar las variables en R hacemos:

## primero codificamos los factores para que estén en el rango -1 y +1 
## (ver # arriba)
coded <- function(x) ifelse(x == x[1], -1, 1)
datos2 <- within(datos2, {
    coded_A = coded(A)
    coded_B = coded(B)
    coded_C = coded(C) 
    coded_D = coded(D)
    
    })
head(datos2, 10)

   RunNumber A B C D RunLabel FiltrationRate coded_D coded_C coded_B coded_A
1          1 − − − −      (1)             45      -1      -1      -1      -1
2          2 + − − −        a             71      -1      -1      -1       1
3          3 − + − −        b             48      -1      -1       1      -1
4          4 + + − −       ab             65      -1      -1       1       1
5          5 − − + −        c             68      -1       1      -1      -1
6          6 + − + −       ac             60      -1       1      -1       1
7          7 − + + −       bc             80      -1       1       1      -1
8          8 + + + −      abc             65      -1       1       1       1
9          9 − − − +        d             43       1      -1      -1      -1
10        10 + − − +       ad            100       1      -1      -1       1

Posteriormente estimamos el modelo de RLM haciendo:

# now fit the regression excluding the interaction term
fitlm <- lm(FiltrationRate ~ coded(A)+coded(C)+coded(D)+coded(A):coded(C)+coded(A):coded(D),datos2)
summary(fitlm)

Call:
lm(formula = FiltrationRate ~ coded(A) + coded(C) + coded(D) + 
    coded(A):coded(C) + coded(A):coded(D), data = datos2)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3750 -1.5000  0.0625  2.9062  5.7500 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         70.062      1.104  63.444 2.30e-14 ***
coded(A)            10.812      1.104   9.791 1.93e-06 ***
coded(C)             4.938      1.104   4.471   0.0012 ** 
coded(D)             7.313      1.104   6.622 5.92e-05 ***
coded(A):coded(C)   -9.063      1.104  -8.206 9.41e-06 ***
coded(A):coded(D)    8.313      1.104   7.527 2.00e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.417 on 10 degrees of freedom
Multiple R-squared:  0.966, Adjusted R-squared:  0.9489 
F-statistic: 56.74 on 5 and 10 DF,  p-value: 5.14e-07

Observe que los para ́metros del modelo ajustado son

\(\hat{\mathbf{\beta}} = (\) 70.0625, 10.8125, 4.9375, 7.3125,-9.0625,8.3125\()\)

que corresponden a las cantidades anteriormente calculadas a partir de los efectos principales.

Validación de supuestos

En el DOE \(2^k\) deben validarse los mismos supuestos sobre el error. En R dicha validación se realiza sobre los residuales del modelo ajustado, que en nuestro corresponde al objeto fitlm.

Para ello utilizaremos la función validar_supuestos() sobre el objeto fitlm. Primero cargamos la función utilizando source()

## cargar la función
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')

Esta función fue creada por el profesor Jorge Vélez, PhD 
<https://jorgeivanvelez.netlify.app/> para la validación 
de los supuestos básicos en RLM. 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a <jvelezv@uninorte.edu.co> 
 
Última modificación:  Agosto 13, 2024 
 

Finalmente aplicamos la función:

## validación de supuestos
validar_supuestos(fitlm)

      normalidad homocedasticidad    independencia 
        "cumple"         "cumple"         "cumple" 

library(ggplot2)

residuals<-residuals(fitlm)
ggplot(data.frame(residuals), aes(sample = residuals)) +
  stat_qq() +
  stat_qq_line(col = "red") +
  labs(title = "Normal Q-Q Plot", x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()

efectos <- c(abs(summary(fitlm[[4]][1])),abs(summary(fitlm[[4]][2])),abs(summary(fitlm[[4]][3])),abs(summary(fitlm[[4]][4])),abs(summary(fitlm[[4]][5])))
ggplot(data.frame(efectos), aes(sample = abs(efectos))) +
  stat_qq(distribution = stats::qnorm) +
  stat_qq_line(distribution = stats::qnorm, col = "blue") +
  labs(title = "Half-Normal Probability Plot", x = "Theoretical Quantiles", y = "Absolute Effects") +
  theme_minimal()

Warning: Removed 24 rows containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 24 rows containing non-finite outside the scale range
(`stat_qq_line()`).

Conclusión: Los resultados indican que los errores del modelo ajustado son independientes, siguen una distribución Normal y tienen varianza constante.

La conclusión global de este análisis es que los factores A, C, D y A:C,A:D tienen un efecto importante sobre la variable FiltrationRate y como los residuales del modelo ajustado cumplen con todos los supuestos, las conclusiones que se deriven de él son válidas para todo el proceso.

Superficie de respuesta

## availability of knitr
if(!require(plotly)) install.packages('plotly')
require(plotly)

# Definir los rangos de factores
A_values <- seq(-1, 1, length.out = 10)
C_values <- seq(-1, 1, length.out = 10)
D_values <- seq(-1, 1, length.out = 10)

# Crear una cuadrícula de valores para A y C
grid <- expand.grid(A = A_values, C = C_values,D=D_values)

# Calcular el término de interacción A:C
grid$AC <- grid$A * grid$C
grid$AD <- grid$A * grid$D

# Calcular respuesta para cada combinación de A y C
grid$FiltrationRate <- summary(fitlm)[[4]][1] + summary(fitlm)[[4]][2] * grid$A + summary(fitlm)[[4]][3]*grid$C +summary(fitlm)[[4]][4]*grid$D + summary(fitlm)[[4]][5]*grid$AC+summary(fitlm)[[4]][5]*grid$AD

# Convertir los valores en una matriz para la visualización
FiltrationRate_matrix <- matrix(grid$FiltrationRate, nrow = 10, ncol = 10)

Warning in matrix(grid$FiltrationRate, nrow = 10, ncol = 10): data length
differs from size of matrix: [1000 != 10 x 10]

fig <- plot_ly(z = ~FiltrationRate_matrix)
fig <- fig %>% add_surface()

fig

Problema 6.1

6.1 In a \(2^4\) factorial design,the number of degrees of freedom for the model,assuming the complete factorial model,is (a)7 (b)5 (c) 6 (d)11 (e) 12 (f) none of the above?

En \(2^k\) el número de grados de libertad totales está dado por la formula \((2^k*n)-1\) donde k es el número de factores \(k=4\) y \(n\) el número de réplicas. Como no se tiene el número de réplicas. La respuesta es la \(f) ninguna de las anteriores\). Si se asumiera una sola réplica el valors de los grados de libertá serian \(glT=2^4-1=15\). Al no estár entre las opciones de respuesta. Se ratificad la respuesta f.

Problema 6.5

An engineer is interested in the effects of cutting speed (A), tool geometry(B), and cutting angle(C)on the life(in hours)of a machine tool.Two levels of each factor are chosen, and three replicates of a \(2^3\) factorial design are run.The results are as follows:

(a)

Estimate the effects of the factors. What effects appear to be large?

Lectura de datos

## read the data set
datos3 <- read.table('C:/Users/FBA/OneDrive - Acesco/0. Personales/2. Profesional/20240416 - Maestria en Analítica de Datos/20240726 Diseño de Experimentos/Homework 2/Datos Problem 6.5.txt', header = TRUE) 

datos3$A<-as.factor(datos3$A)
datos3$B<-as.factor(datos3$B)
datos3$C<-as.factor(datos3$C)

## see first 5 rows of the data
head(datos3, 10)

   Tratamiento A B C lifehours
1          (1) - - -        22
2            a + - -        32
3            b - + -        35
4           ab + + -        55
5            c - - +        44
6           ac + - +        40
7           bc - + +        60
8          abc + + +        39
9          (1) - - -        31
10           a + - -        43

## primero codificamos los factores para que estén en el rango -1 y +1 
## (ver # arriba)
coded <- function(x) ifelse(x == x[1], -1, 1)
datos3 <- within(datos3, {
    coded_A = coded(A)
    coded_B = coded(B)
    coded_C = coded(C) 
    })
head(datos3, 24)

   Tratamiento A B C lifehours coded_C coded_B coded_A
1          (1) - - -        22      -1      -1      -1
2            a + - -        32      -1      -1       1
3            b - + -        35      -1       1      -1
4           ab + + -        55      -1       1       1
5            c - - +        44       1      -1      -1
6           ac + - +        40       1      -1       1
7           bc - + +        60       1       1      -1
8          abc + + +        39       1       1       1
9          (1) - - -        31      -1      -1      -1
10           a + - -        43      -1      -1       1
11           b - + -        34      -1       1      -1
12          ab + + -        47      -1       1       1
13           c - - +        45       1      -1      -1
14          ac + - +        37       1      -1       1
15          bc - + +        50       1       1      -1
16         abc + + +        41       1       1       1
17         (1) - - -        25      -1      -1      -1
18           a + - -        29      -1      -1       1
19           b - + -        50      -1       1      -1
20          ab + + -        46      -1       1       1
21           c - - +        38       1      -1      -1
22          ac + - +        36       1      -1       1
23          bc - + +        54       1       1      -1
24         abc + + +        47       1       1       1

Posteriormente estimamos el modelo de RLM haciendo:

# now fit the regression excluding the interaction term
fitlm <- lm(lifehours~ coded_A*coded_B*coded_C, data=datos3)
summary(fitlm)

Call:
lm(formula = lifehours ~ coded_A * coded_B * coded_C, data = datos3)

Residuals:
   Min     1Q Median     3Q    Max 
-5.667 -3.500 -1.167  3.167 10.333 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              40.8333     1.1211  36.421  < 2e-16 ***
coded_A                   0.1667     1.1211   0.149 0.883680    
coded_B                   5.6667     1.1211   5.054 0.000117 ***
coded_C                   3.4167     1.1211   3.048 0.007679 ** 
coded_A:coded_B          -0.8333     1.1211  -0.743 0.468078    
coded_A:coded_C          -4.4167     1.1211  -3.939 0.001172 ** 
coded_B:coded_C          -1.4167     1.1211  -1.264 0.224475    
coded_A:coded_B:coded_C  -1.0833     1.1211  -0.966 0.348282    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.492 on 16 degrees of freedom
Multiple R-squared:  0.7696,    Adjusted R-squared:  0.6689 
F-statistic: 7.637 on 7 and 16 DF,  p-value: 0.0003977

Observe que los efectos del modelo ajustado son

\(\hat{\mathbf{\beta}} = (\) 81.6666667, 0.3333333, 11.3333333, 6.8333333,-1.6666667,-8.8333333,-2.8333333,-2.1666667\()\)

que corresponden al doble de las cantidades anteriormente calculadas a partir de los efectos principales.El efecto del factor B parece ser más representativo sobre el tiempo de vida de la maquina.

(b)

Use the analysis of variance to confirm your conclusions for part(a).

Construcción de la tabla ANOVA

Inicialmente consideraremos la tabla ANOVA con interacción

## full ANOVA model
fit1 <- aov(lifehours ~ A * B * C, datos3)
summary(fit1)

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1    0.7     0.7   0.022 0.883680    
B            1  770.7   770.7  25.547 0.000117 ***
C            1  280.2   280.2   9.287 0.007679 ** 
A:B          1   16.7    16.7   0.552 0.468078    
A:C          1  468.2   468.2  15.519 0.001172 ** 
B:C          1   48.2    48.2   1.597 0.224475    
A:B:C        1   28.2    28.2   0.934 0.348282    
Residuals   16  482.7    30.2                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se seleccionan con un nivel de significancia \(\alpha = 0.05\) los factores A, B y Cy la interacción A:C. El factor A se incluye por jerarquía en el modelo. Ya que no es estadisticamente significativo, sin emabargo su interacción con C si lo és.

(c)

(c)Write down a regression model for predicting tool life (inhours)based on the results of this experiment.

Construcción del modelo de regresión

# now fit the regression excluding the interaction term
fitlm3 <- lm(lifehours~ coded_A*coded_B*coded_C, data=datos3)
summary(fitlm3)

Call:
lm(formula = lifehours ~ coded_A * coded_B * coded_C, data = datos3)

Residuals:
   Min     1Q Median     3Q    Max 
-5.667 -3.500 -1.167  3.167 10.333 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              40.8333     1.1211  36.421  < 2e-16 ***
coded_A                   0.1667     1.1211   0.149 0.883680    
coded_B                   5.6667     1.1211   5.054 0.000117 ***
coded_C                   3.4167     1.1211   3.048 0.007679 ** 
coded_A:coded_B          -0.8333     1.1211  -0.743 0.468078    
coded_A:coded_C          -4.4167     1.1211  -3.939 0.001172 ** 
coded_B:coded_C          -1.4167     1.1211  -1.264 0.224475    
coded_A:coded_B:coded_C  -1.0833     1.1211  -0.966 0.348282    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.492 on 16 degrees of freedom
Multiple R-squared:  0.7696,    Adjusted R-squared:  0.6689 
F-statistic: 7.637 on 7 and 16 DF,  p-value: 0.0003977

\(\hat{\mathbf{y}}=\) 40.8333333+ 0.1666667\(*x_A\)+ 5.6666667\(*x_B\)+ 3.4166667\(*x_C\)+-4.4166667\(*x_A\)\(x_C\)

(d)

(d)Analyze the residuals.Are there any obvious problems?

Validación de supuestos

## load function
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')

Esta función fue creada por el profesor Jorge Vélez, PhD 
<https://jorgeivanvelez.netlify.app/> para la validación 
de los supuestos básicos en RLM. 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a <jvelezv@uninorte.edu.co> 
 
Última modificación:  Agosto 13, 2024 
 

Ahora validamos el modelo:

## validación de supuestos
validar_supuestos(fitlm3)

      normalidad homocedasticidad    independencia 
     "no cumple"         "cumple"         "cumple" 

los residuos del modelo parecen no cumplir con normalidad, para un nivel de significancia de \(\alpha=0.05\).

## cálculo residuales
r <- residuals(fitlm3)
qqnorm(r)
qqline(r)

plot(fitted(fitlm3), residuals(fitlm3))
abline(h=0)

A partir del objeto r podemos aplicar diferentes funciones para validar los supuestos de Normalidad, Independencia y varianza constante:

## valores p para la validación de supuestos
library(lmtest)

Cargando paquete requerido: zoo

Adjuntando el paquete: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

validacion <- c('Normalidad' = shapiro.test(r)$p.value,
  'Independencia (A)' = dwtest(r ~ A, data = datos3)$p.value,
  'Independencia (B)' = dwtest(r ~ B, data = datos3)$p.value,
  'Independencia (C)' = dwtest(r ~ C, data = datos3)$p.value,
  'Varianza Constante (A)' = bptest(r ~ A, data = datos3)$p.value,
  'Varianza Constante (B)' = bptest(r ~ B, data = datos3)$p.value,
  'Varianza Constante (C)' = bptest(r ~ C, data = datos3)$p.value)
validacion <- round(validacion, 3)
as.list(validacion)

$Normalidad
[1] 0.062

$`Independencia (A)`
[1] 0.769

$`Independencia (B)`
[1] 0.685

$`Independencia (C)`
[1] 0.648

$`Varianza Constante (A).BP`
[1] 0.484

$`Varianza Constante (B).BP`
[1] 0.359

$`Varianza Constante (C).BP`
[1] 0.041

Los residuos, con un nivel de significancias de \(\alpha=0.05\) no parece haber evidencia estadistica que indique que no sigan una distribución normal, o que no sean independiente, o no sean homocedasticos. Exceptuando el factor Ccuyo P_value<0.05

(e)

Análisis de la superficie de respuesta

(e)On the basis of ananalysis of main effect and interaction plots, what coded factor levels of A, B,and C would you recommend using?

# interaction plot

with(datos3, interaction.plot(A, C, lifehours, las = 1))

with(datos3, plot(A,lifehours, las = 1, xlab="A"))

with(datos3, plot(B, lifehours, las = 1, xlab="B"))

with(datos3, plot(C, lifehours, las = 1, xlab="C"))

# Definir los rangos de factores
A_values <- seq(-1, 1, length.out = 10)
C_values <- seq(-1, 1, length.out = 10)
B_values <- seq(-1, 1, length.out = 10)

# Crear una cuadrícula de valores para A y C
grid <- expand.grid(A = A_values, C = C_values,B=B_values)

# Calcular el término de interacción A:C
grid$AC <- grid$A * grid$C

# Calcular respuesta para cada combinación de A y C
grid$lifehours <- summary(fitlm3)[[4]][1] + summary(fitlm3)[[4]][2] * grid$A + summary(fitlm3)[[4]][3] * grid$B +summary(fitlm3)[[4]][4] * grid$C + summary(fitlm3)[[4]][6] * grid$AC

# Convertir los valores en una matriz para la visualización
lifehours_matrix <- matrix(grid$lifehours, nrow = 10, ncol = 10)

Warning in matrix(grid$lifehours, nrow = 10, ncol = 10): data length differs
from size of matrix: [1000 != 10 x 10]

# Crear la visualización de superficie
library(plotly)

fig1 <- plot_ly(z = ~lifehours_matrix, x = ~A_values, y = ~C_values)
fig1 <- fig1 %>% add_surface()
fig1 <- fig1 %>% layout(scene = list(
  xaxis = list(title = "A"),
  yaxis = list(title = "C"),
  zaxis = list(title = "Lifehours")
))
fig1

#maximo de Lifehours
grid$lifehours[which.max(grid$lifehours)]

[1] 54.16667

#Valores máximos
grid$A[which.max(grid$lifehours)]

[1] -1

grid$B[which.max(grid$lifehours)]

[1] 1

grid$C[which.max(grid$lifehours)]

[1] 1

Por lo tanto se sugiere usar valors mínimos de A =-1, máximos de B = +1 y Máximos de C = +1, para obtener un valor máximo del valor de la vida de 54.1666667

Problema 6.15

## read the data set
datos4 <- read.table('C:/Users/FBA/OneDrive - Acesco/0. Personales/2. Profesional/20240416 - Maestria en Analítica de Datos/20240726 Diseño de Experimentos/Homework 2/Datos Problem 6.15.txt', header = TRUE) 

datos4$BottleType<-as.factor(datos4$BottleType)
datos4$Worker<-as.factor(datos4$Worker)

## see first 5 rows of the data
head(datos4, 15)

   BottleType Worker Effort
1       Glass      1     39
2       Glass      1     58
3       Glass      1     44
4       Glass      1     42
5       Glass      2     45
6       Glass      2     35
7       Glass      2     35
8       Glass      2     21
9     Plastic      1     20
10    Plastic      1     16
11    Plastic      1     13
12    Plastic      1     16
13    Plastic      2     13
14    Plastic      2     11
15    Plastic      2     10

#Número de experimentos
NROW(datos4)

[1] 16

## número de réplicas
with(datos4, tapply(Effort, list(BottleType, Worker), length))

        1 2
Glass   4 4
Plastic 4 4

Construcción tabla anova

fit3<-aov(Effort~BottleType*Worker,data=datos4)
summary(fit3)

                  Df Sum Sq Mean Sq F value   Pr(>F)    
BottleType         1 2626.6  2626.6  57.912 6.25e-06 ***
Worker             1  248.1   248.1   5.469   0.0375 *  
BottleType:Worker  1   60.1    60.1   1.324   0.2722    
Residuals         12  544.3    45.4                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los efectos de la interacción entre los factores no es significativa. Por lo cual se plantea el siguiente modelo:

## primero codificamos los factores para que estén en el rango -1 y +1 
## (ver # arriba)
coded <- function(x) ifelse(x == x[1], -1, 1)
datos <- within(datos4, {
    coded_BottleType = coded(BottleType)
    coded_Worker = coded(Worker) 
 
    })
head(datos, 10)

   BottleType Worker Effort coded_Worker coded_BottleType
1       Glass      1     39           -1               -1
2       Glass      1     58           -1               -1
3       Glass      1     44           -1               -1
4       Glass      1     42           -1               -1
5       Glass      2     45            1               -1
6       Glass      2     35            1               -1
7       Glass      2     35            1               -1
8       Glass      2     21            1               -1
9     Plastic      1     20           -1                1
10    Plastic      1     16           -1                1

#Recodificación de variables

fitlm4<-lm(Effort~BottleType+Worker,data=datos4)
summary(fitlm4)

Call:
lm(formula = Effort ~ BottleType + Worker, data = datos4)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.938  -2.188  -0.625   2.031  14.188 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         43.813      2.952  14.840 1.58e-09 ***
BottleTypePlastic  -25.625      3.409  -7.517 4.39e-06 ***
Worker2             -7.875      3.409  -2.310   0.0379 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.818 on 13 degrees of freedom
Multiple R-squared:  0.8263,    Adjusted R-squared:  0.7996 
F-statistic: 30.92 on 2 and 13 DF,  p-value: 1.145e-05

BottleType_values <- Worker_values <-seq(-1,1,length=10)

grid<-expand.grid(BottleType=BottleType_values, Worker=Worker_values)

# Calcular efuerzo para cada combinación de A y C
grid$Effort <- summary(fitlm4)[[4]][1] + summary(fitlm4)[[4]][2] * grid$BottleType + summary(fitlm4)[[4]][3] * grid$Worker

# Convertir los valores de la respuesta en una matriz para la visualización
Effort_matrix <- matrix(grid$Effort, nrow = 10, ncol = 10)

# Crear la visualización de superficie
library(plotly)

fig4 <- plot_ly(z = ~Effort_matrix, x = ~BottleType_values, y = ~Worker_values)
fig4 <- fig4 %>% add_surface()
fig4 <- fig4 %>% layout(scene = list(
  xaxis = list(title = "BottleType"),
  yaxis = list(title = "Worker"),
  zaxis = list(title = "Effort")
))
# Mostrar la gráfica
fig4

Intervalo de confianza

confint(fitlm4)

                      2.5 %      97.5 %
(Intercept)        37.43445  50.1905455
BottleTypePlastic -32.98973 -18.2602674
Worker2           -15.23973  -0.5102674

Los resultados del análisis de varianza y el de los intervalos de confianza coinciden. Ambos intervalos para los factores elegidos, son significativos ya que ninguno de los intervalos incluyen el 0. Con un nivel de significancia de \(\alpha=0.05\) lo cual quiere decir que existe unca combinación para el tipo de de botella y trabajador en el cual los efectos son significativamente diferentes a la media global.