# Cargar librerías necesarias
library(dplyr)
# Importar y guardar archivos planos de textos
library(readr)
# verificación de robots.txt para web scraping
library(polite)
# para manipular datos
library(tidyverse)
# permite realizar web scraping
library(rvest)
## need to load the rsm package
if(!require(rsm)) install.packages('rsm')
require(rsm)
## disponibilidad de lmtest
if(!require(lmtest)) install.packages('lmtest')
require(lmtest)
## disponibilidad de lme4 para componentes de varianza
if(!require(lme4)) install.packages('lme4')
require(lme4)Diseño y Análisis de Experimentos
Resultados Examen #1
CONTEXTO
The Tire Company (TTC) fabrica compuestos para llantas de F1. En la actualidad, TTC desarrolla un nuevo compuesto para la temporada 2025, con miras a superar a su más desafiante competidor.
En el proceso de fabricación existen 3 compuestos \(x_1, x_2 y x_3\) que conforman cada llanta. De acuerdo con las regulaciones, una llanta de F1 debe pesar no más de 13 kilogramos sin incluir el rin y los sensores. Además, se sabe que el costo por cada gramo de material \(j\) , en dólares, es \(c_j = \sqrt{j}\), \(j=1,2,3\)
TTC está convencida de que modificando la composición de la llanta es posible alcanzar un mejor rendimiento (medido en kilómetros recorridos), de tal forma que no sea necesario realizar tantas paradas en pits durante una competencia regular.
Actualmente TTC se encuentra diseñando experimentos para dar respuesta a la problemática. El equipo sabe que la formulación actual permite recorrer \(26 \pm 7\) kilómetros con un set de llantas, y que \(x_1={5, 7, 12}\) \(x_2={3, 12, 20}\) y \(x_3={1, 3, 5, 6}\),todos medidos en libras. Por lo tanto, el resto del peso de la llanta debe completarse con caucho sin procesar.
Punto 1
(10 puntos) Esquematice cómo tomaría los datos necesarios para apoyar la toma de decisiones de TTC. Explique, de manera sencilla, cómo organizaría los datos en un archivo .txt y qué significa una de las filas. Cuál es la unidad experimental?.
# Definir los vectores X1, X2 y X3
X1 <- c(5, 7, 12)
X2 <- c(3, 12, 20)
X3 <- c(1, 3, 5, 6)
# Calcular el número total de combinaciones
num_combinations <- length(X1) * length(X2) * length(X3)
# Inicializar una lista para almacenar los resultados
results <- list()
# Calcular los valores de X4, Peso y Costo_Total para cada combinación
index <- 1
i <-1
for (x1 in X1) {
for (x2 in X2) {
for (x3 in X3) {
# Calcular X4
x4 <- 26 - x1 - x2 - x3
if (x4 < 0){x4 <- 0
i <- 0}
# Calcular Peso
peso <- x1 + x2 + x3
# Calcular Costo_Total
costo_total <- round((sqrt(1) * x1 * 453.592) + (sqrt(2) * x2 * 453.592) + (sqrt(3) * x3 * 453.592),1)
# Agregar los resultados a la lista
results[[index]] <- data.frame(ID = index, X1 = x1, X2 = x2, X3 = x3, X4 = x4, Peso = peso, I=i, Costo_Total = costo_total, km = NA)
index <- index + 1
i <- 1
}
}
}
# Convertir la lista en un data frame
datos <- do.call(rbind, results)
datos$X1 <- as.factor(datos$X1)
datos$X2 <- as.factor(datos$X2)
datos$X3 <- as.factor(datos$X3)
datos$X4 <- as.factor(datos$X4)
# Imprimir el data frame resultante
print(datos) ID X1 X2 X3 X4 Peso I Costo_Total km
1 1 5 3 1 17 9 1 4978.0 NA
2 2 5 3 3 15 11 1 6549.3 NA
3 3 5 3 5 13 13 1 8120.6 NA
4 4 5 3 6 12 14 1 8906.3 NA
5 5 5 12 1 8 18 1 10751.3 NA
6 6 5 12 3 6 20 1 12322.6 NA
7 7 5 12 5 4 22 1 13893.9 NA
8 8 5 12 6 3 23 1 14679.5 NA
9 9 5 20 1 0 26 1 15883.1 NA
10 10 5 20 3 0 28 0 17454.4 NA
11 11 5 20 5 0 30 0 19025.7 NA
12 12 5 20 6 0 31 0 19811.3 NA
13 13 7 3 1 15 11 1 5885.2 NA
14 14 7 3 3 13 13 1 7456.5 NA
15 15 7 3 5 11 15 1 9027.8 NA
16 16 7 3 6 10 16 1 9813.4 NA
17 17 7 12 1 6 20 1 11658.5 NA
18 18 7 12 3 4 22 1 13229.8 NA
19 19 7 12 5 2 24 1 14801.1 NA
20 20 7 12 6 1 25 1 15586.7 NA
21 21 7 20 1 0 28 0 16790.3 NA
22 22 7 20 3 0 30 0 18361.6 NA
23 23 7 20 5 0 32 0 19932.9 NA
24 24 7 20 6 0 33 0 20718.5 NA
25 25 12 3 1 10 16 1 8153.2 NA
26 26 12 3 3 8 18 1 9724.5 NA
27 27 12 3 5 6 20 1 11295.8 NA
28 28 12 3 6 5 21 1 12081.4 NA
29 29 12 12 1 1 25 1 13926.5 NA
30 30 12 12 3 0 27 0 15497.7 NA
31 31 12 12 5 0 29 0 17069.0 NA
32 32 12 12 6 0 30 0 17854.7 NA
33 33 12 20 1 0 33 0 19058.3 NA
34 34 12 20 3 0 35 0 20629.6 NA
35 35 12 20 5 0 37 0 22200.8 NA
36 36 12 20 6 0 38 0 22986.5 NARta.
El esquema sugerido para llevar a cabo el desarrollo del nuevo compuesto es un diseño de experimento factorial fraccionado. Donde \(X_1,X_2,X_3\) como factores con a niveles \(a_1=3,a_2=3,a_4=4\) respectivamente.
Un factor confundente \(X_4\) representa la cantidad de caucho sin procesar. Este factor puede ser medido, pero no controlada, pues dependerá de la combinación de los factores \(X_1,X_2,X_3\) de forma que \(X_4~(X_1,X_2,X_3)=26-X_1-X_2-X_3\)\(X_4 \in [0,26]\)
La variable Peso está determinada por \(Peso~(X_1,X_2,X_3)=X_1+X_2+X_3\)
La variable I está calculado según las restricciones de peso reglamentarias que debe tener una llanta en la F1. \(\forall Peso_i >26 \rightarrow I=0\) y \(\forall Peso_i <26 \rightarrow I=1\) esta variable adicionalmente será un criterio para selección de los fraccionarios
El costo total está dado por la formula Costo_Total \(={453.592*\sqrt1*X_1+453.592*\sqrt2*X_2+453.592*\sqrt2*X_2}\)
La variable de respuesta es Kmen la cual se registrarán los \(Y_{ij}\) Kilometros que resulten de la medición de desgaste de la llanta creada a partir con las combicación de factores identificada con \(ID_i\)
De tal forma que la fila del \(ID_1\) significa la llanta creada a partir de 5 lb del componente \(X_1\), 3 lb del componente \(X_2\), y 1 lb del competente \(X_3\), la cual puede ser rellenada hasta con \(X_4=\) 17 lb de Caucho sin procesar. Este neumático tendríaq un peso de 9lb de los factores controlados. Adicionalmente, la variable \(I_i\)= 9 nos indica que una llanta combinada a partir de esa combinación cumple la reglamentaciones de peso en la F1. El costo total de de producir una llanta con un peso de 1 libra con la combinación \(ID=1\) es de 4978. En la columna Kmse debe diligenciar los kilometros que soporte la llanta hasta ser declarada como desgastada.
Esta tabla debe ser diligenciada y repetirse tantas veces como el tamaño de n.
La unidad experimental sobre las cuales haremos las pruebas es una “Llanta, fabricada a partir de la combinación de compuestos”
Punto 2
(15 puntos) Establezca condiciones bajo las cuales, desde el punto de vista de Diseño de Experimentos, las variables \(X_1\),\(X_2\) y \(X_3\) podrían considerarse como factores de efectos fijos y cuándo como factores de efectos aleatorios. Cuál de los dos escenarios considera que tiene más sentido en el contexto del problema? Justifique su respuesta.
Las variables \(X_1\),\(X_2\)y \(X_3\), podrían ser consideradas de efectos fijos o de efectos aleatorios, según: 1. el propósito del estudio. 2. la relación existente entre dichas variables y la variable de respuesta \(y\). 3. La forma en como sean determinados los niveles de cada factor.
En el caso de Las condiciones principales de un diseño de efectos fijos, el propósito del estudio es saber si la cantidad de los componentes con la que se produzca una llanta tiene un efecto positivo o negativo sobre la cantidad promedio de kilómetros que pueda soportar. Adicionalmente, debemos sospechar o saber que existe una relación entre los factores sobre la media de la respuesta. Debemos también tener interés en conocer los efectos de ciertos niveles específicos.
Por otro lado, en las condiciones de un diseño de efecto aleatoríos, el propósito del estudio es conocer si los factores tienen un impacto sobre la varianza de la distribución de la variable de respuesta. Debe sospecharse o saber que los factores elegidos tienen algún resultado sobre la mejora o decrimento de la precisión en los resultados de la variable de respuesta. Otra condición importante es que los niveles para cada factor deben ser elegidos de manera aleatoria de un pull de niveles existentes.
Para efectos de este experimento, el rendimiento actual es de \(26 \pm 7\) kilómetros. el propósito es alcanzar un mejor redimiento. Es decir Aumentar el promedio de kilómetros recorridos, aumentar el valor de 26. Al revisar la literatura y la información disponible sobre el impacto de los componentes de fabricación de las llantas, se puede evidenciar que el rendimiento de los neumáticos dependerá de los compuestos sobre los cuales fue fabricado. También existen unos niveles elegidos en particular a estudiar para cada efecto (no elegidos aleatoriamente). Por lo cual, el escenario de estudiar efectos fijos, tiene más sentido, en el contexto del problema.
Punto 3
(10 puntos) Considere únicamente las variables \(X_1\) y \(X_2\). Diseñe un experimento de dos factores en efectos fijos teniendo en cuenta que (1) en promedio, cada réplica tarda 2 horas en realizarse; (2) sólo disponemos de 240 horas para tomar una decisión sobre cuáles compuestos utilizar; y (3) la potencia mínima requerida es del 80%. Nota: Una réplica corresponde a un set de llantas.
El diseño de experimento propuesto es del tipo: \(3^k\) donde k es el número de factores estudiados. Para este caso \(K=2\)
Por lo tanto el número total de réplicas está dado por \(N=n*a*k\). Donde \(n:\) número de réplicas realizadas sobre la unidad experimental. \(a:\) número de niveles de cada factor y \(k:\) número de factores.
También se debe considerar que un set de llantas tiene 4 unidades experimentales. Se asume también que en un mismo set de llantas, cada una de las llantas están fabricadas con la misma combinación de compenentes y que el rendimiento de cada llanta para un mismo set se pueden medir de manera individual.
La unidad experimental tde referencia es la llanta
La disposición de la recolección de datos será la siguiente:
# Definir los vectores X1, X2
X1 <- c(5, 7, 12)
X2 <- c(3, 12, 20)
# Calcular el número total de combinaciones
num_combinations <- length(X1) * length(X2)
# Inicializar una lista para almacenar los resultados
results <- list()
# Calcular los valores de X4, Peso y Costo_Total para cada combinación
index <- 1
i <-1
for (x1 in X1) {
for (x2 in X2) {
# Calcular X4
x4 <- 26 - x1 - x2
if (x4 < 0){x4 <- 0
i <- 0}
#Calcular Costo_Total del set
costo_total <- 4*round((sqrt(2) * x1 * 453.592) + (sqrt(2) * x2 * 453.592),1)
# Agregar los resultados a la lista
results[[index]] <- data.frame(ID = index, X1 = x1, X2 = x2, I=i, Costo_Set = costo_total, km1 = NA,km2=NA,Km3=NA,km4=NA)
index <- index + 1
i <- 1
}
}
# Convertir la lista en un data frame
datos <- do.call(rbind, results)
datos$X1 <- as.factor(datos$X1)
datos$X2 <- as.factor(datos$X2)
# Imprimir el data frame resultante
print(datos) ID X1 X2 I Costo_Set km1 km2 Km3 km4
1 1 5 3 1 20527.2 NA NA NA NA
2 2 5 12 1 43620.4 NA NA NA NA
3 3 5 20 1 64147.6 NA NA NA NA
4 4 7 3 1 25659.2 NA NA NA NA
5 5 7 12 1 48752.0 NA NA NA NA
6 6 7 20 0 69279.6 NA NA NA NA
7 7 12 3 1 38488.4 NA NA NA NA
8 8 12 12 1 61581.6 NA NA NA NA
9 9 12 20 0 82108.8 NA NA NA NACada ID corresponde a un set de llantas. La variable Inos indica si es reglamentario o no realizar dicha combinación de los niveles entre los factores \(X_1 y X_2\) dependiendo del peso total de la llanta \(\forall Peso_i >26 \rightarrow I=0\) y \(\forall Peso_i <= 26 \rightarrow I=1\)), por lo cual se deberían excluir del diseño del experimento las combinaciones de los \(ID= {6, 9}\). No realizar pruebas sobre esa combinación de niveles.
El costo del set está cálculado según los gramos de cada factor que compenen cada ID y el valor en gramo para cada componente.
\(km_1, km_2, km_3 y km_4\) serán los valores de rendimiento individuales para cada llanta en el set de llantas probado para la réplica de la combicación ID.
Por lo tanto la tabla de colección de datos sería:
#Eliminamos las combinaciones que no cumplen la restricción del peso según reglamento de F1.
datos <-datos[-c(6,9),]
print(datos) ID X1 X2 I Costo_Set km1 km2 Km3 km4
1 1 5 3 1 20527.2 NA NA NA NA
2 2 5 12 1 43620.4 NA NA NA NA
3 3 5 20 1 64147.6 NA NA NA NA
4 4 7 3 1 25659.2 NA NA NA NA
5 5 7 12 1 48752.0 NA NA NA NA
7 7 12 3 1 38488.4 NA NA NA NA
8 8 12 12 1 61581.6 NA NA NA NANecesitamos sugerir un tamaño de efecto. Para esto consultaremos un poco de la información disponible del competidor principal, con el fin de considerár un tamaño del efecto (bajo, medio o alto) según la relación existente entre la composición de una llanta y su rendimiento :
A través de su Pagina Web encontramos información la realción existente entre los tipos de compuestos y el rendimiento.
El neumático C1 del año pasado se ha convertido en el C0 de 2023. Se trata del neumático más duro de la gama, que estará nominado para los circuitos que más energía sacan de los neumáticos. Será el compuesto preferido si el circuito tiene un asfalto muy abrasivo, si las temperaturas son especialmente altas y si existen numerosas zonas a lo largo de la pista donde los esfuerzos a los que se ven sometidos los neumáticos son elevados.compuesto nuevo para 2023 se ubica entre el C1 y el C2 del año pasado. Creado para reducir la diferencia de rendimiento entre los anteriores dos compuestos más duros de nuestra gama, el C1 se basa en el C2, lo que proporciona una menor degradación y al mismo tiempo ofrece un rendimiento sólido. El C1 se considera una alternativa al C0: los dos no serán seleccionados juntos para ninguna carrera esta temporada.tercer compuesto más duro sigue siendo muy adecuado para los circuitos más rápidos, calientes y abrasivos. Los compuestos más duros también se nominan a veces para nuevos circuitos, con el fin de proporcionar una selección conservadora cuando las circunstancias aún son relativamente desconocidas.tres compuestos nominados para cada carrera, este se usará en todos los eventos por su extremada versatilidad, ya que a veces se usará como el compuesto más duro, a veces como el compuesto medio y, a veces, como el compuesto más blando. Con un excelente equilibrio entre rendimiento y versatilidad, se adapta bien a una amplia gama de circunstancias.compuesto está diseñado para funcionar bien en circuitos estrechos y sinuosos, donde se requiere un calentamiento bastante rápido para alcanzar de inmediato el máximo rendimiento. Este es otro neumático que se usa mucho a lo largo de la temporada, lo que hace adaptable el C4 a varias pistas diferentes.los neumáticos más blandos de la gama, diseñados para los circuitos más lentos con bajo desgaste y degradación en los que se requiere el máximo agarre mecánico de la goma. Estos se usan normalmente en circuitos callejeros, sobre todo en Mónaco, y también en circuitos donde el asfalto es excepcionalmente suave.”
El análisis del texto publicado por el competidor en su página web sugiere que existe un efecto grande entre el tipo de componente usado y el rendimiento final de la llanta.
Un tamaño de efecto (f) grande a detectar. según Cohen (1998) sugiere que valores de (f) de 0.1, 0.25 y 0.4 representan tamaños de efecto pequeño, mediano y grande, respectivamente.
Por lo tanto, utilizaremos \(f=0.4\), y un nivel de significancia de \(\alpha=0.05\)
## availability of pwr
if(!require(pwr)) install.packages('pwr')Cargando paquete requerido: pwrrequire(pwr)
#Cálculo de la potencia
pwr.anova.test(k=length(datos[[1]]),f=0.4,sig.level=0.05,power=0.8)
Balanced one-way analysis of variance power calculation
k = 7
n = 13.09214
f = 0.4
sig.level = 0.05
power = 0.8
NOTE: n is number in each groupn <-ceiling(pwr.anova.test(k=length(datos[[1]]),f=0.4,sig.level=0.05,power=0.8)[[2]])Para alcanzar una potencia de \(power=0.8\) es necesario realizar al menos 14 réplicas, lo que daría un total de pruebas necesarias sobre la unidad experimental de \(N=\) 98
Dado que cada set de llanta probado cuenta con 4 unidades experimentales. El total de set de llantas a probar sería de \(N/4=\) 24.5. Si probar cada ser de llantas toma 2 horas, tenemos un tiempo de experimento de \(T_{exp}=\) 49 horas.
Conclusión: debido a que cada réplica correspondere a un set compuesto de 4 llantas, esto nos permite reducir el tiempo para realizar las pruebas. Dándonos como resultado un total de 49 horas para obtener los resultados de las pruebas. El resto del tiempo dispobible se deberá utilizar en el análsis de la información, para posteriormente tomar decisión. Debido a que el tiempo de experimentación es menor a 240 horas podemos decir que es viable realizar la experimentación
#Cálculamos el número máximo de replicas de set de llantas que se pueden probar en el tiempo disponible
Time_max <- 240
Time_replica <- 2
Num_max_set <- Time_max/Time_replicaPunto 4
Por decisión de TTC, se fijó x3 = 5 y emplearon 9 configuraciones de x1 y x2 con n=3 réplicas. Los datos obtenidos se encuentran aquí
## read the data set
datos <- read.table('https://www.dropbox.com/scl/fi/w4kki7slj3fu26ie2s14a/ttc.txt?rlkey=b02uzksyy1sdkpryjp18wixlh&e=1&dl=1',
header = TRUE)
## definir factores
datos$x1 <- as.factor(datos$x1)
datos$x2 <- as.factor(datos$x2)
datos$x3 <- as.factor(datos$x3)
## see first 3 rows of the data
head(datos, 30) id x1 x2 x3 distance rep
1 1 5 3 5 30.525 1
2 2 5 3 5 30.923 2
3 3 5 3 5 30.439 3
4 4 7 3 5 30.787 1
5 5 7 3 5 30.743 2
6 6 7 3 5 31.156 3
7 7 12 3 5 31.784 1
8 8 12 3 5 31.213 2
9 9 12 3 5 31.634 3
10 10 5 12 5 29.501 1
11 11 5 12 5 29.433 2
12 12 5 12 5 29.137 3
13 13 7 12 5 29.570 1
14 14 7 12 5 30.109 2
15 15 7 12 5 29.944 3
16 16 12 12 5 30.572 1
17 17 12 12 5 30.925 2
18 18 12 12 5 30.499 3
19 19 5 20 5 28.674 1
20 20 5 20 5 28.958 2
21 21 5 20 5 28.305 3
22 22 7 20 5 28.674 1
23 23 7 20 5 28.420 2
24 24 7 20 5 28.707 3
25 25 12 20 5 29.520 1
26 26 12 20 5 30.008 2
27 27 12 20 5 29.841 3El experimento corresponde a un diseño con \(3^k\) con \(k=2\) y un bloque \(x_3=5\),
El número total de experimentos es
## how many runs?
NROW(datos)[1] 27y el número de réplicas por cada tratamiento puede obtenerse haciendo:
## replicates per combination
with(datos, tapply(distance, list(x1, x2, x3), length)), , 5
3 12 20
5 3 3 3
7 3 3 3
12 3 3 34a. Análisis gráfico
Gráficamente los resultados del experimento pueden representarse utilizando gráficos de interacción:
## profile plot x1 vs x2
with(datos, interaction.plot(x1 ,x2, distance), fun = function(x) sum(x),
lty = 1, las = 1,
col = 2:4, type = 'x2', lwd = 2,
pch = 16, xlab = 'Factor x1',
ylab = "Totales", ylim = x3(5))
Conclusión: De acuerdo con este gráfico, parece que existe un efecto principal de cada factor sobre la media de la distancia recorrida. Siendo \(x2=3\) el nivel en el que mejor media de distancia se optiene, para cualquier valor de \(x_1\). La interacción entre ambos compuestos parece no ser importante.
4b. Construcción de la Tabla ANOVA
Realizamos la tabla anova para el modelo completo incluyendo interacción para validar entre \(x_1 y x_2\) para validar a través del valor P, con \(\alpha=0.05\), si es significativa o no la interacción. Se debe excluir del anova la variable \(x3\) dado a que al fijarse en un sola nivel, no se considera un factor, ya que no hay variabilidad en él.
## full ANOVA model
fit <- aov(distance ~ x1+x2+x1:x2, data = datos)
summary(fit) Df Sum Sq Mean Sq F value Pr(>F)
x1 2 6.264 3.132 49.911 4.53e-08 ***
x2 2 18.211 9.105 145.102 7.91e-12 ***
x1:x2 4 0.394 0.099 1.571 0.225
Residuals 18 1.130 0.063
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1conclusión: En este caso, no tiene sentido considerar la interacción x1 y x2. Dado que \(P_{value}(x1:x2)>0.05\). Los efectos de los factores x1 y x2, son significativos. De manera que existe al menos un x1 y un x2 en el que la distancia recorrida son diferentes a la media.
Teniendo en cuenta este resultado, recalcularemos la tabla ANOVA, esta vez sin incluir el término de interacción.
## ANOVA model with no interaction and second order
fit_no_interaction <- aov(distance ~ x1+x2, data = datos)
summary(fit_no_interaction) Df Sum Sq Mean Sq F value Pr(>F)
x1 2 6.264 3.132 45.22 1.61e-08 ***
x2 2 18.211 9.105 131.45 5.82e-13 ***
Residuals 22 1.524 0.069
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Superficie de respuesta
Una vez detectamos que los factores de interés tienen un efecto sobre la variable respuesta \(distance\), es importante cuantificar dicho efecto. Para ello usamos un modelo de Regresión Lineal Múltiple (RLM).
Por simplicidad, nos enfocarmeos en la boquilla tipo 1, es decir, que sólo utilizaremos los datos del experimento correspondientes x1 x1 = 1. Adicionalmente, recodificaremos las variables utilizando la función recode``EnR` primero codificamos las variables
## availability of car package
if(!require(car)) install.packages('car')
require(car)
## recoding
datos$x1 <- as.numeric(as.character(recode(datos$x1, "5 = 5; 7 = 7; 12 = 12")))
datos$x2 <- as.numeric(as.character(recode(datos$x2, "3 = 3; 12 = 12; 20 = 20") ))y luego estimamos el modelo, se realiza evaluación del modelo lineal, incluyendo factores de segundo orden, con el fin de validar su significancia dentro del modelo.
# now fit the regression excluding the double interaction term
type_2 <- lm(distance ~ x1+x2,
data = subset(datos,x3==5))
summary(type_2)
Call:
lm(formula = distance ~ x1 + x2, data = subset(datos, x3 == 5))
Residuals:
Min 1Q Median 3Q Max
-0.46479 -0.14581 -0.04142 0.16140 0.43310
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 30.074176 0.164025 183.351 < 2e-16 ***
x1 0.163201 0.016643 9.806 7.20e-10 ***
x2 -0.118264 0.007056 -16.762 9.45e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2546 on 24 degrees of freedom
Multiple R-squared: 0.9402, Adjusted R-squared: 0.9352
F-statistic: 188.6 on 2 and 24 DF, p-value: 2.106e-15Al realizar el modelo lineal optenemos un \(R^2_{adj}=\) 0.940165
# now fit the regression excluding the triple interaction term
type_1 <- lm(distance ~ x1 + x2 + I(x1^2),
data = subset(datos,x3==5))
summary(type_1)
Call:
lm(formula = distance ~ x1 + x2 + I(x1^2), data = subset(datos,
x3 == 5))
Residuals:
Min 1Q Median 3Q Max
-0.4782 -0.1302 -0.0548 0.1769 0.3996
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 30.569302 0.737863 41.430 < 2e-16 ***
x1 0.033589 0.188963 0.178 0.860
x2 -0.118264 0.007134 -16.577 2.77e-14 ***
I(x1^2) 0.007456 0.010826 0.689 0.498
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2574 on 23 degrees of freedom
Multiple R-squared: 0.9414, Adjusted R-squared: 0.9337
F-statistic: 123.1 on 3 and 23 DF, p-value: 2.616e-14Al realizar el modelo lineal incluyendo terminos de orden cuadrático optenemos ajuste mejor con \(R^2_{adj}=\) 0.940165. Por lo tanto usaremos los siguientes parámetros para los cálculos posteriores:
\(\mathbf{\hat{\beta}} = (\) 30.5693017 , 0.0335889, -0.118264, 0.0074556\()\)
x1 partir del modelo de RLM ajustado, contenido en el objeto type_1, 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ínimo de \(distance\) diferentes valores de \(x_1\) y \(x_2\). Recordemos que \(x_j\in \{-1, 0, 1\}\), \(j=1,2,3\).
# grid-based optimization
h <- 1000
x1grid <- seq(1, 50, length = h)
x2grid <- seq(1, 50, length = h)
grid <- expand.grid(x1 = x1grid, x2 = x2grid)
## now calculate ypred and transform to matrix
y_pred <- predict(type_1, newdata = grid)
y_pred <- matrix(y_pred, ncol = length(x1grid), nrow = length(x1grid))x1 continuación graficamos la superficie de respuesta:
## check availability of the plot3D package
if(!require(plot3D)) install.packages('plot3D')
require(plot3D)
## 3D plot
persp3D(x1grid, x2grid, y_pred,
xlab = "\n Compuesto X1",
ylab = "\n Compuesto X2",
zlab = "\n Distance",
colkey = FALSE,
ticktype = 'detailed',
main = "Compuesto X3 = 5")
Por otro lado, los contornos de la superficie son:
## countourplot
contour2D(y_pred, x1grid, x2grid, xlab = "Compuesto X1",
ylab = "Compuesto X2", colkey = FALSE, las = 1, col = c(2:10),
main = "X3 = 5")
Validación de supuestos
Validamos los supuestos :
## x3álculo de residuales
r <- residuals(type_1)
## valor p
c('p_normalidad' = shapiro.test(r)$p.value,
'p_independencia' = car:::durbinWatsonTest(type_1)$p) p_normalidad p_independencia
0.2130285 0.5100000 Para validar el supuesto de varianza constante podemos hacerlo para cada factor x1 través de la prueba de Bartlett:
## valor p de la prueba de Bartlett por factor
c('p_x1' = bartlett.test(r ~ datos$x1)$p.value,
'p_x2' = bartlett.test(r ~ datos$x2)$p.value) p_x1 p_x2
0.9974153 0.9452106 ## cargando al 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
Ahora validamos los supuestos:
## validación de supuestos x1 partir de un modelo de RLM
validar_supuestos(type_1) normalidad homocedasticidad independencia
"cumple" "cumple" "cumple" Conclusión validación supuestos: El modelo ajustado cumple con todos los supuestos. Por tanto, las conclusiones derivadas de él son válidas para la población de la que provienen los datos.El modelo nos permitirá realizar pronósticos y estimado de los optimos de distancia según los compuestos elegidos.
Adicionalmente, la compración de los modelos sugiere que existe una de segúndo orden con la variable de \(x_1\). Y su inclusión nos asegura un mejor \(R^2_{adj]}\)
4c. Analisis de varianza efecto aleatorios
(15 puntos) Realice el análisis anterior bajo en elfoque de efectos aleatorios y estime las componentes de varianza. Observa alguna diferencia? Concluya.
Construcción de la Tabla ANOVA
Realizamos la tabla anova para las varianzas \(x_1 y x_2\) para validar a través del valor P, con \(\alpha=0.05\), si es significativa o no la interacción. Se debe excluir del anova la variable \(x3\) dado a que al fijarse en un sola nivel, no se considera un factor, ya que no hay variabilidad en él.
## full ANOVA model
fit2 <- aov(distance ~ Error(x1*x2), data = datos)
summary(fit2)
Error: x1
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 6.232 6.232
Error: x2
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 18.21 18.21
Error: x1:x2
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 0.09102 0.09102
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 23 1.465 0.06368 MSS1 <- broom::tidy((aov(distance ~ Error(x1*x2), data = datos)))[[5]][1]
MSS2 <- broom::tidy((aov(distance ~ Error(x1*x2), data = datos)))[[5]][2]
MSS3 <- broom::tidy((aov(distance ~ Error(x1*x2), data = datos)))[[5]][3]
MSS4 <- broom::tidy((aov(distance ~ Error(x1*x2), data = datos)))[[5]][4]A partir de estos resultados, los estadísticos de prueba para cada factor pueden obtenerse como:
\(F_{x1 x x2}=\) \({MS_{x1 x x2}/MSE}=\) 0.09/0.06\(=\) 1.43
\(F_{x1}=\) \({MS_{x1}/MSE}=\) 6.23/0.06\(=\) 97.88
\(F_{x2}=\) \({MS_{TCopperContent}/MSE}=\) 18.21/0.06\(=\) 285.98
\(F_{Error}=\) \({MSE}=\) 0.06
Estimación de las Componentes de Varianza
A continuación se muestra la sintaxis para estimar los componentes de varianza.
## full variance components models
vcm <- lmer(distance ~ 1 + (1|x1) + (1|x2) + (1|x1:x2),
data = datos)
summary(vcm)Linear mixed model fit by REML ['lmerMod']
Formula: distance ~ 1 + (1 | x1) + (1 | x2) + (1 | x1:x2)
Data: datos
REML criterion at convergence: 24.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.6075 -0.5672 -0.1302 0.6171 1.4117
Random effects:
Groups Name Variance Std.Dev.
x1:x2 (Intercept) 0.01194 0.1093
x2 (Intercept) 1.00075 1.0004
x1 (Intercept) 0.33704 0.5806
Residual 0.06275 0.2505
Number of obs: 27, groups: x1:x2, 9; x2, 3; x1, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 30.0000 0.6705 44.74Por lo tanto,
\[ \begin{eqnarray} \hat{\sigma}^2_{\text{x1}} &=& 1.001 \\\nonumber \hat{\sigma}^2_{\text{x2}} &=& 0.337 \\\nonumber \hat{\sigma}^2_{\text{x1 * x2}} &=& 0.012\\\nonumber \hat{\sigma}^2_{\text{Error}} &=& 0.063 \end{eqnarray} \]
Finalmente,
\(\hat{\sigma}^2_y =\) 1.001 + 0.337 + 0.012 + 0.063 = 1.413
De esta variabilidad total, el factor x1 contribuye con un 70.84%, x2 con un 23.85%, la interacción x1:x2 con un 0.85% y finalmente el Error un un 4.46%.
Pruebas de hipótesis
Una forma probar si las componentes de varianza asociadas a x2, x1 y x2:x1 son estadísticamente significativas para una probabilidad de Error Tipo I del 5% es a través del cálculo de intervalos de confianza. Para ello utilizaremos la función confint():
## intervalos de confianza del 95%
confint(vcm, oldNames = FALSE) 2.5 % 97.5 %
sd_(Intercept)|x1:x2 0.0000000 0.4245282
sd_(Intercept)|x2 0.4370429 2.5707700
sd_(Intercept)|x1 0.2522892 1.7594997
sigma 0.1865734 0.3588956
(Intercept) 28.5099878 31.4900898Así las cosas,
\({\sigma^2}_{x1 x x2}\in(\) 0\(^2,\) 0.42 \(^2) \rightarrow {\sigma^2}_{x1 x x2} \in (\) 0 \(,\) 0.18\()\)
\({\sigma^2}_{x2}\in(\) 0.44\(^2,\) 2.57 \(^2) \rightarrow {\sigma^2}_{x2} \in (\) 0.19 \(,\) 6.61 \()\)
\({\sigma^2}_{x1}\in(\) 0.25\(^2,\) 1.76 \(^2) \rightarrow {\sigma^2}_{x1} \in (\) 0.06 \(,\) 3.1 \()\)
\({\sigma^2}_{Error}\in(\) 0.19\(^2,\) 0.36 \(^2) \rightarrow {\sigma^2}_{Error} \in (\) 0.03 \(,\) 0.13 \()\)
A partir de este resultado podemos concluir que ni x2 ni la interacción x1:x2 son fuentes significativas de la variabilidad de y.
Validación de Supuestos
Con el fin de determinar si las conclusiones obtenidas hasta ahorap pueden o no extenderse para toda la población, es necesario validar los supuestos.
Los residuales del modelo puede estimarse haciendo:
## cálculo residuales
r <- residuals(vcm)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
validacion <- c('Normalidad' = shapiro.test(r)$p.value,
'Independencia (x2)' = dwtest(r ~ x2, data = datos)$p.value,
'Independencia (x1)' = dwtest(r ~ x1, data = datos)$p.value,
'Varianza Constante (x2)' = bptest(r ~ x2, data = datos)$p.value,
'Varianza Constante (x1)' = bptest(r ~ x1, data = datos)$p.value)
validacion <- round(validacion, 3)
as.list(validacion)$Normalidad
[1] 0.232
$`Independencia (x2)`
[1] 0.678
$`Independencia (x1)`
[1] 0.719
$`Varianza Constante (x2).BP`
[1] 0.904
$`Varianza Constante (x1).BP`
[1] 0.966Conclusión: Como los valores \(p\) de todas las pruebas son \(>0.05\), podemos afirmar que los resultados obtenidos son generalizables para la población de la que proviene la muestra.
Punto 5
(30 puntos) Cuáles son los valores de \(x_1\) y \(x_2\) tal que la distancia recorrida esté en el rango \(32 \pm 3\) Km? Calcule el costo de un set de llantas para estas condiciones.
Utilizando la superficie de respuesta anterior limitamos los valores para los compuestos \(X_1 y X_2\) de forma que las isocuantas estén entre 29 y 35 Kilómetros. Adicionalmente, incluimos la restricción reglamentaria de forma que \(X_1+X_2+X_3 \leq 26\) como \(X_3=5 \rightarrow X_1+X_2 \leq 21\)
De forma que:
Teniendo en cuenta la restricción \(X_1+X_2 \leq 21\)
Podemos sustituir \(X_2= 21-X_1\) en \(Y\). Obteniendo: Para encontrar el valor máximo posible de y, debemos maximizar la función:
y = 30.569302 + 0.033589x1 - 0.118264x2 + 0.007456x1² + 0.2574
Sujeta a las restricciones:
\(x_1 + x_2 ≤ 21\), \(x_1>0\) \(x_2>0\)
Evaluamos la funciones y restricciones para \(x_1=0\)
\((0) + x_2 = 21\)
y = 30.569302 + 0.033589*(0) - 0.118264(21)+ 0.007456(0)² + 0.2574 y = 34.82
Evaluamos las funciones y restricciones para \(x_2=0\)
\(x_1 + (0) = 21\)
y = 30.569302 + 0.033589*(21) - 0.118264(0)+ 0.007456(21)² + 0.2574 y = 28.3432
Estos son los máximos y mínimos que podemos tener para nuestro modelo de regresión, sujeto a las restricciones existentes. Lo primero que podemos evidenciar es que no es posible encontrar una combinación de los componentes \(x_1 y x_2\) tal que la distancia resultante sea de 35. máximo será de \(Y_{max}(0,21)= 34.82\), por otro lado \(Y_{min}(21,0)= 28.3432\) nos da indicio de que es viable encontrar una combinación de los componentes \(x_1 y x_2\) tal que \(Y=29\).
Para encontrar el valor de x1 y x2 cuando y = 29, podemos sustituir y en la ecuación original:
29 = 30.569302 + 0.033589x1 - 0.118264x2 + 0.007456x1² + 0.2574
Como x1 + x2 ≤ 21, podemos reescribir x2 como:
x2 =21-x1
Sustituimos x2 en la ecuación:
21 -x1 = 0.0630454 x1^2 + 0.284017 x1 + 15.446
Simplificamos la ecuación:
0 = 0.0630454 x1^2+1.284017-5.554
Resolvemos la ecuación cuadrática para x1:
x1 ≈ 3.665 y -24.03 (no viable)
Encontramos valor de x2 cuando x1 = 3.66
x2=21-(3.66) = 17.334
# grid-based optimization
h <- 1000
x1grid <- seq(0, 23, length = h)
x2grid <- seq(0, 17.334, length = h)
grid <- expand.grid(x1 = x1grid, x2 = x2grid)
y_pred <- predict(type_1, newdata = grid)
y_pred <- matrix(y_pred, ncol = length(x1grid), nrow = length(x1grid))Graficamos la nueva superficie de posibles respuestas limitando \(X_1 \in [0,22.5]\) y \(X_1 \in [0,13.5]\).
## countourplot
contour2D(y_pred, x1grid, x2grid, xlab = "Compuesto X1",
ylab = "Compuesto X2", colkey = FALSE, las = 1,
main = "X3 = 5")
curva <-function(x1grid) 21-x1grid
curve(curva,from =0,to=26,add = TRUE)
legend("topright", legend = "X2+X1 <= 21", lty = 1, col = "black")
# Sombrea la zona que cumple la restricción de pesaje.
x1 <- seq(0, 26, by = 0.1)
y1 <- curva(x1)
polygon(c(x1, rev(x1)), c(y1, rep(0, length(x1))), col = rgb(0,0,0,alpha=0.1), border = NA)
Conclusión: el costo del set de llantas oscilará según la combinación de componentes que se elijan dentro de las isocuantas existentes dentro del rango habilitado para x1 y x2. \(Y_{min}(3.66,17.334)=29\) y \(Y_{max}(0,21)=34.82\)
Por lo tanto:
x1<-c(3.66,0)
x2<-c(17.334,21)
x3<-c(5,5)
costo_llanta <- round((sqrt(1) * x1* 453.592) + (sqrt(2) * x2 * 453.592) + (sqrt(3) * x3 * 453.592),1)
costo_set<-c(4*costo_llanta[1],4*costo_llanta[2])
costo_set[1] 66830.8 69596.8El costo total de un set de llantas (4 llantas) ocsilará entre 6.68308^{4} y 6.95968^{4}, que no superen un peso máximo de 26 lb, y que tengan un rendimiento entre 29 y 35 km.