library(ggplot2)
library(plotly)
#Datos
observacion <- c(1:25)
tamano <- c(80,30,50,90,70,60,120,80,100,50,40,70,90,20,110,100,30,50,90,110,30,90,40,80,70)
horas<- c(399,121,221,376,361,224,546,352,353,157,160,252,389,113,435,420,212,268,377,421,273,468,244,342,323)
#Modelo
modelo <- lm(horas~tamano)
#Ajustados y residuales
ajustados <- modelo$fitted.values
residuales <- modelo$residuals
#Dataframe
datos <- data.frame(observacion, ajustados, residuales)
# Para líneas horizontales en el gráfico
maximo <-max(residuales)
minimo <- min(residuales)
#Grafico
grafico <- ggplot(data=datos, aes(x=ajustados, y=residuales)) + geom_point()
grafico <- grafico + xlab("Ajustados") + ylab("Residuales") + ggtitle("Verificación Homocedasticidad")
grafico <- grafico + geom_hline (yintercept = c(maximo+5,minimo-5) )
grafico <- grafico + geom_hline(yintercept=0,linetype="dashed", color="blue")
grafico <- grafico + theme_classic()
ggplotly(grafico)Verificación de la adecuación del modelo - Homocedasticidad
En es importante examinar la idoneidad del modelo para los datos antes de usar, de manera práctica, el modelo. En este apartado del curso, se discutirán, métodos gráficos y formales para verificar la adecuación del modelo con respecto a los principales supuestos. Los principales supuestos a tener en cuenta son:
- Los errores están normalmente distribuidos.
- El término del error \(\epsilon\) tiene varianza \(\sigma^2\) constante.
- Los errores no están autocorrelacionados.
- La relación entre \(y\) y los regresores \(x_k\) es lineal, al menos, aproximadamente.
Si se cumple la normalidad y la no autoccorelación, implícitamente los errores son variables aleatorias independientes.
El estudio de los supuestos anteriores se realiza principalmente sobre los residuales del modelo.
Residuales
Los residuales \(e_i\) son definidos como la diferencia entre el valor observado de \(y_i\) y el valor ajustado o estimado \(\hat{y_i}\), como se muestra en la Ecuación 1.
\[\begin{align} e_i = y_i - \hat{y_i} \end{align} \tag{1}\]
Los residuales pueden ser visto como el error observado a diferencia del verdadero error \(\epsilon_i\) desconocido en el modelo de regresión (Ecuación 2).
\[\begin{align} \epsilon_i = y_i - E(y_i) \end{align} \tag{2}\]
Para el modelo de regresión lineal simple, Ecuación 3.
\[\begin{align} y_i = \beta_0 + \beta_1x_i + \epsilon_i \end{align} \tag{3}\]
Los términos del error \(\epsilon\) son asumidos como variables aleatorias normales independientes, con media \(0\) y varianza \(\sigma^2\) constante.
Si el modelo es adecuado para los datos en cuestión, los residuales observados \(e_i\) deberían reflejar las propiedades asumidas para \(\epsilon\);. Esta es la idea básica que subyace en el análisis residual, un medio muy útil para examinar la idoneidad de un modelo estadístico.
Propiedades de los residuales
Media
La media de los \(n\) residuales \(e_i\) para el modelo de regresión lineal simple, Ecuación 3, está dada en la Ecuación 4:
\[\begin{align} \bar{e} = \frac{\sum_{i=1}^n e_i}{n} =0 \end{align} \tag{4}\]
Como la media de los residuales \(e_i\) es siempre 0, no presenta ninguna información sobre si los errores verdaderos \(\epsilon_i\) tienen valor esperado $E[_i] =0 $
Varianza
La varianza para los \(n\) residuales \(e_i\) para el modelo de regresión lineal simple, Ecuación 3, es definida como sigue (Ecuación 5):
\[\begin{align} s^2 &= \frac{\sum_{i=n}^n (e_i - \bar{e})^2}{n-2} \\ \\ s^2 &= \frac{\sum_{i=n}^n e_i^2}{n-2} \\ \\ s^2 &= \frac{SSE}{n-2} \\ \\ s^2 &= MSE \end{align} \tag{5}\]
Residuales semiestudentizados
En ocasiones, resulta útil estandarizar los residuales para el análisis. Dado que la desviación típica de los términos de error \(\epsilon_i\); es \(\sigma\), que se estima por \(\sqrt{MSE}\), es natural considerar la siguiente forma de estandarización:
\[\begin{align} e_i^* &= \frac{e_i - \bar{e}}{\sqrt{MSE}} \\ \\ e_i^* &= \frac{e_i}{\sqrt{MSE}} \end{align} \tag{6}\]
\(e_i^*\) de la Ecuación 7 son llamados residuales estudentizados o semiestudentizados.
Verificación gráfica del supuesto de homocedasticidad.
Una forma de verificar el supuesto de varianza constante (los términos del error tienen la misma varianza) es realizando un gráfico de dispersión con las siguientes variables:
- Eje X \(\rightarrow\) ajustados o predichos \(\hat{y_{i}}\)
- Eje Y \(\rightarrow\) residuales \(e_{i}\)
Si los puntos de este gráfico de dispersión se distribuyen de manera aleatoria en una banda horizontal (sin ningún patrón claro y contundente), entonces es señal de que se cumple el supuesto de que los tratamientos tienen igual varianza.
Verificación gráfica ejemplo Toluca Company
La empresa Toluca Company fabrica equipos de refrigeración, así como muchas piezas de repuesto. En el pasado, una de las piezas de recambio se fabricaba periódicamente en lotes de distintos tamaños. Cuando se emprendió un programa de mejora de costos, los responsables de la empresa quisieron determinar el tamaño de lote óptimo para producir esta pieza. La producción de esta pieza implica la puesta a punto del proceso de producción (que debe realizarse sea cual sea el tamaño del lote) y operaciones de mecanizado y montaje. Un dato clave para que el modelo determinara el tamaño de lote óptimo fue la relación entre el tamaño de lote y las horas de trabajo necesarias para producir el lote. Para determinar esta relación, se utilizaron datos sobre el tamaño de lote y las horas de trabajo de 25 series de producción recientes. Las condiciones de producción fueron estables durante el periodo de seis meses en el que se realizaron las 25 series y se esperaba que siguieran siendo las mismas durante los tres años siguientes, correspondiente al periodo de planificación para el que se estaba llevando a cabo el programa de mejora de costos.
Los datos sobre tamaño de lote y las horas de trabajo de encuentran en la Tabla 1
| Observación \(i\) | \(x_i\): Tamaño de lote | \(y_i\): Horas de trabajo |
|---|---|---|
| 1 | 80 | 399 |
| 2 | 30 | 121 |
| 3 | 50 | 221 |
| 4 | 90 | 376 |
| 5 | 70 | 361 |
| 6 | 60 | 224 |
| 7 | 120 | 546 |
| 8 | 80 | 352 |
| 9 | 100 | 353 |
| 10 | 50 | 157 |
| 11 | 40 | 160 |
| 12 | 70 | 252 |
| 13 | 90 | 389 |
| 14 | 20 | 113 |
| 15 | 110 | 435 |
| 16 | 100 | 420 |
| 17 | 30 | 212 |
| 18 | 50 | 268 |
| 19 | 90 | 377 |
| 20 | 110 | 421 |
| 21 | 30 | 273 |
| 22 | 90 | 468 |
| 23 | 40 | 244 |
| 24 | 80 | 342 |
| 25 | 70 | 323 |
Se realiza el cálculo de los ajustados \(\hat{y_i}\) y residuales \(e_i\), Tabla 2:
| Observación \(i\) | \(\hat{y_i}\): Ajustados | \(e_i\): Residuales |
|---|---|---|
| 1 | 347.9820 | 51.0179798 |
| 2 | 169.4719 | -48.4719192 |
| 3 | 240.8760 | -19.8759596 |
| 4 | 383.6840 | -7.6840404 |
| 5 | 312.2800 | 48.7200000 |
| 6 | 276.5780 | -52.5779798 |
| 7 | 490.7901 | 55.2098990 |
| 8 | 347.9820 | 4.0179798 |
| 9 | 419.3861 | -66.3860606 |
| 10 | 240.8760 | -83.8759596 |
| 11 | 205.1739 | -45.1739394 |
| 12 | 312.2800 | -60.2800000 |
| 13 | 383.6840 | 5.3159596 |
| 14 | 133.7699 | -20.7698990 |
| 15 | 455.0881 | -20.0880808 |
| 16 | 419.3861 | 0.6139394 |
| 17 | 169.4719 | 42.5280808 |
| 18 | 240.8760 | 27.1240404 |
| 19 | 383.6840 | -6.6840404 |
| 20 | 455.0881 | -34.0880808 |
| 21 | 169.4719 | 103.5280808 |
| 22 | 383.6840 | 84.3159596 |
| 23 | 205.1739 | 38.8260606 |
| 24 | 347.9820 | -5.9820202 |
| 25 | 312.2800 | 10.7200000 |
Se realiza gráfico de dispersión \(\hat{y_i}\) vs \(e_i\), Figura 1
Según el Figura 1 no se debería sospechar sobre la violación del supuesto de homocedasticidad.
Verificación formal - analítica supuesto de homocedasticidad mediante Test Breuch Pagan.
El test de Breuch - Pagan para verificar la varianza constante asume que los términos del error son independientes y se distribuyen normalmente y que la varianza del término del error \(\epsilon_i\) denotada como \(\sigma^2\), está relacionada con el nivel de de \(x\), de la siguiente manera, (Ecuación 7):
\[\begin{align} log_e \sigma_i^2 = \gamma_0 + \gamma_1x_i \end{align} \tag{7}\]
1. Planteamiento de la hipótesis.
Note que Ecuación 7 implica que \(\sigma_i^2\) incrementa o disminuye con el cambio en \(x\), dependiendo del signo de \(\gamma_1\). La varianza constante se da en el caso que \(\gamma_1 = 0\). Por lo anterior se pueden plantear las hipótesis de la siguiente forma (Ecuación 8):
\[\begin{align} H_0: \gamma_1 &= 0 \\ H_1: \gamma_1 &\ne 0 \end{align} \tag{8}\]
2. Estadístico de prueba
Existen dos formas de cálculo del estadístico de prueba, el estadístico de prueba \(BP\) propuesto por los autores del test. y el estadístico de prueba \(X_{BP}^2\) propuesto por Cook y Weisberg como mejora del test.
Forma 1. Estadístico de prueba Breuch Pagan \(BP\).
El estadístico de prueba \(BP\) se calcula de la siguiente manera, Ecuación 9:
\[\begin{align} BP = nR^{*2} \end{align} \tag{9}\]
Donde:
- \(R^{*2}\) corresponde al coeficiente de determinación de la regresión siendo \(e_i^2\) la variable respuesta y \(x\) la variable regresora.
Forma 2. Estadístico de prueba modificado Cook and Weisberg \(X_{BP}^2\).
El estadístico de prueba \(X_{BP}^2\) se calcula de la siguiente manera, Ecuación 10:
\[\begin{align} X_{BP}^2 =\frac{\frac{SSR^*}{2}}{\left[ \frac{SSE}{n} \right]^2} \end{align} \tag{10}\]
Donde:
- \(SSR^*\) corresponde a la suma de cuadrados de la regresión siendo \(e_i^2\) la variable respuesta y \(x_i\) la variable regresora.
- \(SSE\) corresponde a la suma de cuadrados del error siendo \(y_i\) la variable respuesta y \(x_i\) la variable regresora.
3. Estadístico teórico
Los estadísticos de prueba siguen una distribución Chi-Cuadrado con \(1\) grado de libertad (en general son \(k\) grados de libertad iguales a la cantidad de variables regresoras)(Ecuación 11 - Ecuación 12), esto es:
\[\begin{align} BP\sim X^2_{\alpha,~1} \end{align} \tag{11}\]
\[\begin{align} X_{BP}^2\sim X^2_{\alpha,~1} \end{align} \tag{12}\]
Rechazo \(H_0\) si:
\[\begin{align} X_{BP}^2 > X^2_{\alpha,~1} \end{align}\]
Verificación formal en R (indirecta) con Test Breuch-Pagan ejemplo Toluca Company usando estadístico de prueba \(X_{BP}^2\).
La empresa Toluca Company fabrica equipos de refrigeración, así como muchas piezas de repuesto. En el pasado, una de las piezas de recambio se fabricaba periódicamente en lotes de distintos tamaños. Cuando se emprendió un programa de mejora de costos, los responsables de la empresa quisieron determinar el tamaño de lote óptimo para producir esta pieza. La producción de esta pieza implica la puesta a punto del proceso de producción (que debe realizarse sea cual sea el tamaño del lote) y operaciones de mecanizado y montaje. Un dato clave para que el modelo determinara el tamaño de lote óptimo fue la relación entre el tamaño de lote y las horas de trabajo necesarias para producir el lote. Para determinar esta relación, se utilizaron datos sobre el tamaño de lote y las horas de trabajo de 25 series de producción recientes. Las condiciones de producción fueron estables durante el periodo de seis meses en el que se realizaron las 25 series y se esperaba que siguieran siendo las mismas durante los tres años siguientes, correspondiente al periodo de planificación para el que se estaba llevando a cabo el programa de mejora de costos.
Los datos sobre tamaño de lote y las horas de trabajo de encuentran en la . Se realizará división de grupos según tamaño del lote.
#Datos
tamano <- c(80,30,50,90,70,60,120,80,100,50,40,70,90,20,110,100,30,50,90,110,30,90,40,80,70)
horas<- c(399,121,221,376,361,224,546,352,353,157,160,252,389,113,435,420,212,268,377,421,273,468,244,342,323)
#Cálculo SSE para modelo y~x
modelo1 <- lm(horas~tamano)
anova1 <- aov(modelo1)
summary(anova1) Df Sum Sq Mean Sq F value Pr(>F)
tamano 1 252378 252378 105.9 4.45e-10 ***
Residuals 23 54825 2384
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SSE<-54825
#Cálculo SSE para modelo e^2~x
residuales <- anova1$residuals
residuales_2 <-residuales^2
modelo2 <- lm(residuales_2~tamano)
anova2<- aov(modelo2)
summary(anova2) Df Sum Sq Mean Sq F value Pr(>F)
tamano 1 7896142 7896142 1.091 0.307
Residuals 23 166395896 7234604 SSR <- 7896142
#Estadístico de Prueba
x_bp <- (SSR/2)/((SSE/25)^2)
x_bp[1] 0.820933#Estadístico teórico
qchisq(1-0.05,1)[1] 3.841459De los resultados en R se obtiene que, para \(\alpha = 0.05\):
\[\begin{align} X_{BP}^2 = 0.820933 \ngtr t_{0.95,~1} = 3.841459 \end{align}\]
Por lo que no existe evidencia suficiente para rechazar \(H_0\), concluyendo que en el modelo \(log_e \sigma_i^2 = \gamma_0 + \gamma_1x_i\) el parámetro \(\gamma_1 = 0\), la varianza no depende de los \(x_i\), por lo tanto, se cumple el supuesto de homocedasticidad.
Verificación formal en R (indirecta) con Test Breuch-Pagan ejemplo Toluca Company usando estadístico de prueba \(BP\).
#Datos
tamano <- c(80,30,50,90,70,60,120,80,100,50,40,70,90,20,110,100,30,50,90,110,30,90,40,80,70)
horas<- c(399,121,221,376,361,224,546,352,353,157,160,252,389,113,435,420,212,268,377,421,273,468,244,342,323)
#Modelo y~x
modelo1 <- lm(horas~tamano)
#Modelo e^2~x
residuales <- modelo1$residuals
residuales_2 <-residuales^2
modelo2 <- lm(residuales_2~tamano)
# Coeficiente de determinación modelo 2
R2 <- summary(modelo2)$r.squared
#Estadístico de Prueba
BP <- length(tamano)*R2
BP[1] 1.132602#Estadístico teórico
qchisq(1-0.05,1)[1] 3.841459Verificación formal en R (directa) con Test Breuch-Pagan ejemplo Toluca Company para estadístico de prueba BP.
#Datos
tamano <- c(80,30,50,90,70,60,120,80,100,50,40,70,90,20,110,100,30,50,90,110,30,90,40,80,70)
horas<- c(399,121,221,376,361,224,546,352,353,157,160,252,389,113,435,420,212,268,377,421,273,468,244,342,323)
# Modelos de datos
modelo<- lm(horas~tamano)
#Breuch - Pagan Test
library(lmtest)
bptest(modelo)
studentized Breusch-Pagan test
data: modelo
BP = 1.1326, df = 1, p-value = 0.2872#Estadístico teórico
qchisq(1-0.05,1)[1] 3.841459