2. Regresión lineal simple (en espanol)

2. Regresión lineal simple¶

La Toluca Company fabrica equipos de refrigeración y las piezas de repuesto que los mantienen funcionando. Una pieza se había producido durante mucho tiempo en lotes, o tandas, del tamaño que en cada momento pareciera conveniente. Cuando la empresa inició una revisión para reducir costos, sus ingenieros quisieron resolver una pregunta básica: ¿qué tamaño de lote es el más barato de producir? Cada corrida de producción conlleva un costo fijo de preparación que no depende del tamaño del lote, más el trabajo de maquinado y ensamble que crece con el número de piezas. Para encontrar el punto óptimo, los ingenieros necesitaban primero un ingrediente: ¿cómo dependen las horas de trabajo que consume una corrida del tamaño del lote?

Extrajeron los registros de 25 corridas recientes realizadas bajo condiciones estables a lo largo de un periodo de seis meses. Cada registro es un par: el tamaño del lote (unidades producidas) y las horas de trabajo que tomó la corrida. Los tamaños de lote son todos múltiplos de diez, una conveniencia de programación, y van de 20 a 120 unidades. La dispersión de estos 25 puntos es la imagen inicial de todo el capítulo.

Dos cosas resaltan en la Figure 1. Los lotes más grandes toman más horas, y la nube de puntos está lo bastante cerca de una recta como para que una recta sea un resumen razonable. Pero ¿cuál recta? Un ingeniero podría poner una regla sobre el gráfico y trazar una a ojo, y un segundo ingeniero trazaría una un poco distinta. Este capítulo reemplaza la regla por una regla en el otro sentido: una receta precisa y defendible para la única mejor recta, un recuento honesto de cuánto confiar en ella, y los supuestos que hacen válida toda la empresa. Al terminar podrás ajustar esa recta, leer su pendiente en las unidades del problema (horas por unidad de tamaño de lote), decir qué tan incierta es la pendiente, y demostrar que, bajo las condiciones establecidas, ningún método competidor lo hace mejor.

El flujo de trabajo del modelado¶

Cada capítulo de este libro es una vuelta por el mismo ciclo de cinco etapas. Nombrar las etapas ahora te da un mapa para llevar a lo largo de todo el curso, y cada capítulo posterior comienza diciendo a qué etapa sirve.

• PREGUNTAR. Fija la pregunta y los datos que pueden responderla. Aquí: ¿cómo dependen las horas de trabajo del tamaño del lote, y qué nos pueden decir 25 corridas de producción?

• EXPLORAR. Grafica y resume antes de ajustar nada. La dispersión en la Figure 1 es esta etapa; es lo que nos dice que vale la pena siquiera intentar una recta.

• AJUSTAR. Estima el modelo. La mayor parte de este capítulo vive aquí, convirtiendo la recta de regla y ojo en la única recta de mínimos cuadrados.

• VERIFICAR. Pregunta si el modelo ajustado es confiable, usando los residuos y los supuestos que hay detrás de ellos. Empezamos esto en la Sección 2.3 y dedicamos todo el Capítulo 9 a ello.

• USAR. Interpreta las estimaciones, predice y decide. Leer la pendiente en horas por unidad de tamaño de lote es el primer sabor; la inferencia y la predicción llenan el Capítulo 3.

Como muestra la Figure 2, las etapas forman un ciclo, no una calle de un solo sentido: una VERIFICACIÓN fallida te devuelve a EXPLORAR o AJUSTAR con un mejor modelo. La regresión lineal simple es el viaje completo más pequeño alrededor de este ciclo, y por eso el curso empieza aquí.

2.1 El modelo de regresión lineal simple¶

Intuición¶

La regresión parte de una idea sencilla: una cosa tiende a moverse con otra, y queremos precisar esa relación en números. Aquí describimos cómo una respuesta $Y$ (horas de trabajo de Toluca) depende de un solo predictor $X$ (tamaño del lote). Ningún modelo sensato afirma que las horas son una función exacta del tamaño del lote: dos corridas del mismo tamaño tomaron cantidades de tiempo distintas, por el desgaste de las herramientas, el turno de trabajadores, el clima y cientos de pequeñas cosas que nadie registró. Así que dividimos cada observación en dos piezas: una parte sistemática que el tamaño del lote explica, y una parte aleatoria que no.

La parte sistemática es una recta. En promedio, una corrida de tamaño $X$ toma cierta cantidad base de tiempo más un número fijo de horas extra por cada unidad adicional. La parte aleatoria es una perturbación que empuja las horas reales por encima o por debajo de ese promedio. Suponemos que la perturbación no tiene una tendencia incorporada a ser positiva o negativa, tiene el mismo tamaño típico sin importar el tamaño del lote, y no se coordina de una corrida a la siguiente.

Fórmula¶

El modelo de regresión lineal simple (simple linear regression model) (Definición 2.1) para las observaciones $i = 1, \dots, n$ es

• $Y_i$ es la respuesta observada en el $i$-ésimo caso (horas de trabajo de la corrida $i$).

• $X_i$ es el predictor en el $i$-ésimo caso (tamaño del lote de la corrida $i$), tratado como una constante conocida y fija. Algunos libros lo llaman variable independiente o variable explicativa; este libro dice predictor en todo el texto.

• $\beta_0$ es el parámetro de intercepto: la respuesta media cuando $X = 0$.

• $\beta_1$ es el parámetro de pendiente: el cambio en la respuesta media por un aumento de una unidad en $X$.

• $\varepsilon_i$ es el error aleatorio en el $i$-ésimo caso, la cantidad por la que $Y_i$ se aparta de su media.

Los parámetros $\beta_0$ y $\beta_1$ son números fijos que no conocemos y queremos estimar. El modelo impone cuatro condiciones sobre los errores:

En palabras: cada error promedia a cero (la recta es correcta en promedio), todo error tiene la misma varianza $\sigma^2$ (dispersión constante, llamada homocedasticidad), y errores distintos no están correlacionados (la suerte de una corrida no dice nada sobre la de otra). Nota lo que todavía no se supone: ninguna forma particular para la distribución del error. La normalidad entra solo en 2.6 Máxima verosimilitud bajo errores normales, y solo cuando queremos construir intervalos y pruebas.

Como $X_i$ es una constante fija y $\beta_0 + \beta_1 X_i$ es por lo tanto también constante, las condiciones sobre los errores se transfieren directamente a las respuestas:

La cantidad $E\{Y_i\}$ es la respuesta media (mean response) al nivel $X_i$. La función $E\{Y\} = \beta_0 + \beta_1 X$ que da la respuesta media en cualquier $X$ es la función de regresión (regression function) (Definición 2.2). Es una recta, y estimarla es todo el juego.

La Figure 3 es la imagen mental que hay que llevar a lo largo del capítulo: una recta media, con una distribución de posibles valores de $Y$ apilados en cada $X$, cada pila del mismo ancho.

2.2 Mínimos cuadrados desde los primeros principios¶

Intuición¶

Tenemos una nube de puntos y queremos la única recta que mejor se ajusta. “Mejor” necesita una definición. Los mínimos cuadrados la definen penalizando los fallos: para una recta candidata, mide la brecha vertical de cada punto a la recta, eleva al cuadrado las brechas para que pasarse y quedarse corto cuenten ambos como errores y los fallos grandes duelan de forma desproporcionada, y súmalas. La mejor recta es la que hace este total, la suma de los errores al cuadrado, lo más pequeño posible. La Figure 4 muestra la idea con tres rectas candidatas a través de los datos de Toluca.

Fórmula¶

Escribe una recta candidata con intercepto $b_0$ y pendiente $b_1$. Su suma de errores al cuadrado (sum of squared errors) (Definición 2.3) es

• $Y_i - b_0 - b_1 X_i$ es la brecha vertical del punto $(X_i, Y_i)$ a la recta candidata.

• $Q$ es una función de los dos números $b_0, b_1$ que podemos elegir; los datos están fijos.

Las estimaciones de mínimos cuadrados $b_0$ y $b_1$ son los valores que minimizan $Q$. Para escribirlas de forma compacta, define las tres sumas básicas de cuadrados y productos cruzados alrededor de las medias:

donde $\bar{X}$ y $\bar{Y}$ son las medias muestrales del predictor y la respuesta. Las dos derivaciones siguientes establecen los minimizadores, enunciados aquí primero.

Derivación por cálculo¶

Demostración (por cálculo). $Q(b_0, b_1)$ es una suma de cuadrados, así que es suave y está acotada inferiormente por cero. En un mínimo, ambas derivadas parciales se anulan. Deriva:

Igualando cada una a cero y dividiendo entre -2 se obtienen las ecuaciones normales (normal equations):

Desarrolla las sumas. La primera ecuación se convierte en $\sum Y_i = n b_0 + b_1 \sum X_i$, y dividiendo entre $n$ da $\bar{Y} = b_0 + b_1 \bar{X}$, de modo que

Sustituye esto en la segunda ecuación normal $\sum X_i Y_i = b_0 \sum X_i + b_1 \sum X_i^2$. Reemplazando $b_0$ y agrupando los términos de $b_1$,

El lado izquierdo es igual a $\sum X_i Y_i - n \bar{X}\bar{Y} = \sum (X_i - \bar X)(Y_i - \bar Y) = S_{xy}$, y el corchete del lado derecho es igual a $\sum X_i^2 - n\bar{X}^2 = \sum (X_i - \bar X)^2 = S_{xx}$. Por lo tanto

Este es el único punto estacionario, y es un mínimo por el criterio de las segundas derivadas para una función de dos variables (la misma prueba que usarías en cálculo para una variable, extendida a dos). Las tres segundas derivadas parciales de $Q$ son constantes: $Q_{b_0 b_0} = 2n$, $Q_{b_1 b_1} = 2\sum X_i^2$ y $Q_{b_0 b_1} = 2\sum X_i$. Forma la cantidad de prueba $D = Q_{b_0 b_0}\,Q_{b_1 b_1} - (Q_{b_0 b_1})^2 = 4\big(n\sum X_i^2 - (\sum X_i)^2\big) = 4n\,S_{xx}$. Como $D > 0$ siempre que los $X_i$ no sean todos iguales, y $Q_{b_0 b_0} = 2n > 0$, el punto estacionario es un mínimo local. Como $Q$ es una suma de cuadrados (así que está acotada inferiormente por cero) y este es su único punto estacionario, ese mínimo local es también el mínimo global. (El Capítulo 7 revisita este mismo hecho en forma matricial, una vez que tengas las herramientas del Capítulo 6.) $\blacksquare$

Derivación por una identidad algebraica¶

El argumento de cálculo encuentra el minimizador pero deja la geometría un poco oculta. Aquí hay una segunda demostración que no usa derivadas y muestra directamente que ninguna recta le gana a la recta de mínimos cuadrados. Descansa sobre dos hechos que demostraremos en 2.3 Valores ajustados y las propiedades de los residuos: escribiendo los valores ajustados como $\hat{Y}_i = b_0 + b_1 X_i$ y los residuos como $e_i = Y_i - \hat{Y}_i$, los residuos cumplen $\sum e_i = 0$ y $\sum X_i e_i = 0$.

Demostración (por una identidad algebraica). Toma cualquier recta competidora con intercepto $a$ y pendiente $c$. Divide la brecha de cada punto a esa recta pasándola por el valor ajustado de mínimos cuadrados:

Eleva al cuadrado y suma. La suma de errores al cuadrado de la recta competidora es

El término cruzado se anula:

Así que para toda recta competidora,

con igualdad solo cuando $\sum g_i^2 = 0$, es decir cuando $g_i = 0$ para todo $i$, lo cual (dado que los $X_i$ no son todos iguales) fuerza $a = b_0$ y $c = b_1$. La recta de mínimos cuadrados es el minimizador único. $\blacksquare$

Las dos ecuaciones normales tienen una imagen limpia. Cada una es una recta en el plano $(b_0, b_1)$, y las estimaciones de mínimos cuadrados están donde se cruzan, en el fondo del tazón de SSE. La Figure 5 muestra el tazón y la Figure 6 lo muestra desde arriba.

R¶

Ajustar el modelo en R usa lm, que lee una fórmula response ~ predictor. Primero lee los datos y míralos.

toluca <- read.csv("data/toluca.csv")
dim(toluca)
head(toluca, 3)

[1] 25  2
  lotsize hours
1      80   399
2      30   121
3      50   221

Poder reproducir lm a partir de las sumas crudas $S_{xx}$ y $S_{xy}$ importa: significa que puedes calcular una regresión en cualquier lugar, verificar la salida del software y, más adelante, extender la idea cuando no haya un botón lm para el modelo que necesitas.

Vale la pena calcular los conjuntos de datos pequeños a mano una vez, para que las fórmulas dejen de ser cajas negras. El archivo toluca_mini.csv es un recorte de seis filas de los datos de Toluca, lo bastante pequeño para hacerlo en papel. La Figure 7 muestra sus puntos y sus desviaciones respecto a las medias, la materia prima de $S_{xx}$ y $S_{xy}$.

Antes de aceptar las fórmulas por confianza, dedica un minuto a intentar ganarles a mano.

Qué observar. Ningún par de posiciones lleva $Q$ por debajo de $54{,}825.5$, y las que más se acercan marcan 62.5 y 3.57, las estimaciones de mínimos cuadrados. Prueba esto. Deja la pendiente en 3.57 y mueve solo el intercepto: la suma de residuos cruza cero justo en la posición donde $Q$ toca fondo, que es la primera ecuación normal ocurriendo frente a ti. Volver a la Sección 2.2.

2.3 Valores ajustados y las propiedades de los residuos¶

Intuición¶

Una vez fijada la recta, cada punto observado se divide en dos partes. El valor ajustado (fitted value) $\hat{Y}_i$ (Definición 2.5) es donde la recta dice que debería estar el punto; el residuo (residual) $e_i = Y_i - \hat{Y}_i$ es qué tan lejos cae el punto real de él. Los residuos son las sobras del modelo, la parte de $Y$ que el tamaño del lote no explicó. La Figure 9 los dibuja como los palitos verticales que conectan cada punto con la recta. No son solo decoración: los mínimos cuadrados fuerzan a los residuos a cumplir identidades algebraicas exactas, y esas identidades son a la vez una herramienta de demostración y una forma práctica de atrapar un error en un ajuste.

Fórmula¶

Los valores ajustados y los residuos son

• $\hat{Y}_i$ es la respuesta media estimada en $X_i$, la altura de la recta ajustada allí.

• $e_i$ es el residuo, la brecha vertical de la $Y_i$ observada hacia abajo hasta la recta.

No confundas $e_i$ con el error del modelo $\varepsilon_i = Y_i - (\beta_0 + \beta_1 X_i)$. El error es la brecha a la recta verdadera y desconocida; el residuo es la brecha a nuestra recta estimada. Nunca vemos $\varepsilon_i$, pero siempre podemos calcular $e_i$, y es nuestro mejor sustituto. La Figure 10 pone las dos brechas lado a lado para que la diferencia sea imposible de pasar por alto.

Derivación (propiedades de los residuos)¶

Demostración. Todas ellas se siguen de las dos ecuaciones normales de 2.2 Mínimos cuadrados desde los primeros principios.

1. Los residuos suman cero. La primera ecuación normal es exactamente $\sum (Y_i - b_0 - b_1 X_i) = \sum e_i = 0$.

2. El predictor es ortogonal a los residuos. La segunda ecuación normal es exactamente $\sum X_i (Y_i - b_0 - b_1 X_i) = \sum X_i e_i = 0$.

3. Los valores ajustados son ortogonales a los residuos. Usando las propiedades 1 y 2, $\sum \hat{Y}_i e_i = \sum (b_0 + b_1 X_i) e_i = b_0 \sum e_i + b_1 \sum X_i e_i = 0$.

4. Los valores ajustados reproducen el total de $Y$. De la propiedad 1, $\sum e_i = \sum (Y_i - \hat{Y}_i) = 0$, así que $\sum \hat{Y}_i = \sum Y_i$, y dividiendo entre $n$, el valor ajustado medio es igual a $\bar{Y}$.

5. La recta pasa por el centro de los datos. Poniendo $X = \bar{X}$ en la recta ajustada y usando $b_0 = \bar{Y} - b_1 \bar{X}$ da $\hat{Y} = b_0 + b_1 \bar{X} = \bar{Y}$. Así que $(\bar{X}, \bar{Y})$ está siempre sobre la recta de mínimos cuadrados. $\blacksquare$

La propiedad 1 es por qué los palitos verticales de la Figure 9 se equilibran. Las propiedades 1 y 2 fueron los ingredientes que hicieron desaparecer el término cruzado en la demostración algebraica de 2.2 Mínimos cuadrados desde los primeros principios, así que las identidades pagan el préstamo que tomamos allí.

R y Python¶

Un gráfico de los residuos contra el predictor es la primera imagen que hay que trazar después de cualquier ajuste. Si se cumplen los supuestos de varianza constante y de recta del modelo, los residuos deberían moverse alrededor de cero sin tendencia y sin abanicarse, como en la Figure 11. Una curva señalaría que la relación no es recta; un embudo que se ensancha señalaría varianza no constante. El Capítulo 9 construye un juego de herramientas de diagnóstico completo sobre este único hábito.

Las identidades del Teorema 2.6 son exactas, así que la forma más rápida de creerlas es intentar romper una.

Qué observar. Por mucho que destroces los datos, $\sum e_i$ y $\sum X_i e_i$ se quedan en cero y la media de $\hat{Y}$ sigue igual a $\bar{Y}$. Prueba esto. Arrastra una corrida hasta lo más alto del gráfico: los residuos se vuelven enormes, el ajuste queda inservible, y las identidades no se mueven, y por eso verifican tu aritmética y nunca tu modelo. Volver a la Sección 2.3.

2.4 Estimar la varianza del error¶

Intuición¶

La pendiente y el intercepto describen la recta, pero el modelo tiene una tercera incógnita: $\sigma^2$, la varianza de los errores, que controla qué tan estrechamente los puntos abrazan la recta. Una $\sigma^2$ pequeña significa que las predicciones pueden ser afiladas; una grande significa que incluso una recta perfecta deja mucha incertidumbre. No podemos ver los errores $\varepsilon_i$, pero los residuos $e_i$ son sus sustitutos visibles, así que estimamos $\sigma^2$ a partir de la dispersión de los residuos. La Figure 13 muestra cómo se ve esa dispersión como una banda alrededor de la recta ajustada.

Fórmula¶

La suma de cuadrados del error (error sum of squares) (Definición 2.7) y el estimador de $\sigma^2$ son

• $\mathrm{SSE}$ es el residuo total al cuadrado, la variación sobrante que la recta no explicó.

• $\mathrm{MSE}$ (cuadrado medio del error), también escrito $s^2$, es $\mathrm{SSE}$ dividido entre sus grados de libertad $n-2$; estima $\sigma^2$.

• $s$ es la desviación estándar estimada de los errores, en las unidades de $Y$ (horas).

El divisor es $n-2$, no $n$ ni $n-1$. Dos grados de libertad se gastan estimando los dos coeficientes $b_0$ y $b_1$ que definen la recta desde la que se miden los residuos. La siguiente derivación muestra que este divisor exacto, y ningún otro, hace que $\mathrm{MSE}$ sea correcto en promedio.

Derivación (por qué el divisor es $n-2$)¶

Demostración. Mostramos que $E\{\mathrm{SSE}\} = (n-2)\sigma^2$, de modo que dividir entre $n-2$ da un estimador que es correcto en promedio. Parte de una identidad ordenada para la SSE. Como $\hat{Y}_i = \bar{Y} + b_1(X_i - \bar{X})$ (usando $b_0 = \bar{Y} - b_1 \bar X$),

así que, desarrollando el cuadrado y usando $S_{xy} = b_1 S_{xx}$,

Toma esperanzas de las dos piezas. Para cualquier variable aleatoria, $E\{W^2\} = \operatorname{Var}\{W\} + (E\{W\})^2$.

La pieza de la pendiente. Mostramos en 2.5 Comportamiento muestral y el teorema de Gauss-Markov que $E\{b_1\} = \beta_1$ y $\operatorname{Var}\{b_1\} = \sigma^2 / S_{xx}$. Por lo tanto

La pieza total. Escribe $S_{yy} = \sum Y_i^2 - n\bar{Y}^2$. Cada $Y_i$ tiene media $\mu_i = \beta_0 + \beta_1 X_i$ y varianza $\sigma^2$, así que $E\{Y_i^2\} = \sigma^2 + \mu_i^2$ y $E\{\sum Y_i^2\} = n\sigma^2 + \sum \mu_i^2$. La media muestral $\bar{Y}$ tiene media $\bar{\mu} = \beta_0 + \beta_1 \bar{X}$ y varianza $\sigma^2/n$, así que $E\{n\bar{Y}^2\} = n(\sigma^2/n + \bar{\mu}^2) = \sigma^2 + n\bar{\mu}^2$. Restando,

El corchete final es $\sum (\mu_i - \bar{\mu})^2 = \sum \big(\beta_1 (X_i - \bar{X})\big)^2 = \beta_1^2 S_{xx}$. Así que $E\{S_{yy}\} = (n-1)\sigma^2 + \beta_1^2 S_{xx}$.

Combinar. Restando la pieza de la pendiente de la pieza total,

Por lo tanto $E\{\mathrm{MSE}\} = E\{\mathrm{SSE}\}/(n-2) = \sigma^2$: el divisor $n-2$ es exactamente lo que hace que $\mathrm{MSE}$ sea un estimador insesgado de $\sigma^2$. $\blacksquare$

R y Python¶

El Teorema 2.8 dice que el divisor $n-2$ es el que acierta en promedio, y el muestreo repetido es la manera de ver cómo se asienta esa afirmación.

Qué observar. Los valores individuales de $\mathrm{MSE}$ se dispersan mucho, la mayoría cae entre unos 1200 y unos 4000, y aun así su promedio camina firme hacia 2383.4. Prueba esto. Presiona Tomar 1000 dos veces y compara las dos medias: el divisor $n-2$ cae en $\sigma^2$ mientras que el divisor $n$ se queda cerca de 8 por ciento por debajo por muchas muestras que agregues, y eso es lo que significa el sesgo. Volver a la Sección 2.4.

2.5 Comportamiento muestral y el teorema de Gauss-Markov¶

Intuición¶

Las estimaciones $b_0$ y $b_1$ se calculan a partir de una muestra particular de 25 corridas. Otras 25 corridas, hechas bajo las mismas condiciones, darían estimaciones un poco distintas. Así que $b_0$ y $b_1$ son en sí mismos aleatorios: tienen distribuciones, llamadas distribuciones muestrales, a lo largo de las muestras repetidas que podríamos haber extraído. La Figure 15 lo hace concreto simulando 60 muestras nuevas en los tamaños de lote de Toluca y ajustando cada una: las rectas ajustadas se dispersan alrededor de la recta media verdadera, y se abren más donde $X$ está lejos de $\bar{X}$.

Dos preguntas organizan esta sección. ¿Son $b_0$ y $b_1$ correctos en promedio (insesgados)? ¿Y cuánto varían de muestra a muestra (sus varianzas)? El teorema de Gauss-Markov entonces entrega una recompensa fuerte: entre todos los estimadores sensatos que uno podría inventar, los mínimos cuadrados tienen la varianza más pequeña.

Fórmula¶

Ambos estimadores son insesgados, y sus varianzas son

• $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$ se encoge cuando los errores son pequeños ($\sigma^2$ bajo) o el predictor está muy disperso ($S_{xx}$ grande). Dispersar tus valores de $X$ compra precisión.

• $\operatorname{Var}\{b_0\}$ carga un término extra $\bar{X}^2/S_{xx}$ porque el intercepto vive lejos en $X = 0$, e inclinar la recta la hace pivotar allí.

La Figure 16 hace visible el primer hecho. Cuando los tamaños de lote están apiñados, un poco de ruido inclina la pendiente descontroladamente; cuando están muy dispersos, el mismo ruido apenas la mueve.

Reemplazar la desconocida $\sigma^2$ por su estimación $\mathrm{MSE}$ da las varianzas estimadas $s^2\{b_1\} = \mathrm{MSE}/S_{xx}$ y $s^2\{b_0\} = \mathrm{MSE}\big(1/n + \bar{X}^2/S_{xx}\big)$, cuyas raíces cuadradas son los errores estándar $s\{b_1\}$ y $s\{b_0\}$ que el software imprime.

Derivación (insesgadez y varianzas)¶

Demostración (media y varianza de $b_1$). La clave es escribir $b_1$ como una suma ponderada de las respuestas. Como $\sum (X_i - \bar{X}) \bar{Y} = 0$,

Los pesos $k_i$ dependen solo de los valores fijos de $X$ y cumplen tres hechos, cada uno un cálculo de una línea:

Insesgadez. Usando $E\{Y_i\} = \beta_0 + \beta_1 X_i$,

Varianza. Los $Y_i$ no están correlacionados con varianza común $\sigma^2$, así que las varianzas de una suma ponderada se suman con pesos al cuadrado:

Demostración (media y varianza de $b_0$). Escribe $b_0 = \bar{Y} - b_1 \bar{X} = \sum m_i Y_i$ con $m_i = \tfrac{1}{n} - \bar{X} k_i$. Entonces $\sum m_i = 1$ y $\sum m_i X_i = \bar X - \bar X = 0$, así que $E\{b_0\} = \beta_0 \sum m_i + \beta_1 \sum m_i X_i = \beta_0$. Para la varianza, usando $\sum k_i = 0$ y $\sum k_i^2 = 1/S_{xx}$,

Derivación (Gauss-Markov: $b_1$ es BLUE)¶

Demostración. Considera cualquier estimador lineal $\hat\beta_1^{\ast} = \sum c_i Y_i$ con pesos fijos $c_i$. Su esperanza es $\sum c_i(\beta_0 + \beta_1 X_i) = \beta_0 \sum c_i + \beta_1 \sum c_i X_i$. Para que esto sea igual a $\beta_1$ para todo posible $\beta_0$ y $\beta_1$, los pesos deben cumplir

Escribe $c_i = k_i + d_i$, dividiendo cada peso en el peso de mínimos cuadrados $k_i$ y una desviación $d_i$. Como tanto $c_i$ como $k_i$ cumplen las dos restricciones de insesgadez, las desviaciones cumplen $\sum d_i = 0$ y $\sum d_i X_i = 0$. Ahora el término cruzado desaparece:

Por lo tanto

Como $\sigma^2 \sum d_i^2 \ge 0$, todo competidor lineal insesgado tiene varianza al menos la de $b_1$, con igualdad solo cuando todos los $d_i = 0$, es decir cuando el competidor es $b_1$. La pendiente de mínimos cuadrados es BLUE. El mismo argumento con los pesos $m_i$ muestra que $b_0$ es BLUE. $\blacksquare$

R y Python¶

Las fórmulas de varianza no son solo álgebra; describen el haz que se abre de rectas en la Figure 15. Podemos confirmar tanto la insesgadez como la fórmula de la varianza de la pendiente por simulación, extrayendo 5000 muestras del modelo estimado y reajustando cada una. La Figure 17 muestra el histograma de las 5000 estimaciones de pendiente.

set.seed(4210)
beta0 <- 62.366; beta1 <- 3.5702; sigma <- 48.82
xfix <- toluca$lotsize
Sxx_fix <- sum((xfix - mean(xfix))^2)
b1_sim <- replicate(5000, {
  ysim <- beta0 + beta1 * xfix + rnorm(length(xfix), 0, sigma)
  coef(lm(ysim ~ xfix))[["xfix"]]
})
round(c(mean_of_b1_hats = mean(b1_sim),
        true_beta1 = beta1,
        sd_of_b1_hats = sd(b1_sim),
        formula_se = sqrt(sigma^2 / Sxx_fix)), 4)

mean_of_b1_hats      true_beta1   sd_of_b1_hats      formula_se
         3.5699          3.5702          0.3479          0.3469

rng = np.random.default_rng(4210)
beta0, beta1, sigma = 62.366, 3.5702, 48.82
xfix = toluca["lotsize"].to_numpy()
Sxx_fix = np.sum((xfix - xfix.mean()) ** 2)
b1_hats = np.empty(5000)
for i in range(5000):
    ysim = beta0 + beta1 * xfix + rng.normal(0, sigma, size=xfix.size)
    b1_hats[i] = np.polyfit(xfix, ysim, 1)[0]
print(round(b1_hats.mean(), 4), beta1,
      round(b1_hats.std(ddof=1), 4),
      round(np.sqrt(sigma ** 2 / Sxx_fix), 4))

3.5675 3.5702 0.352 0.3469

El promedio de las 5000 estimaciones de pendiente queda justo sobre $\beta_1 = 3.5702$ (insesgadez), y su desviación estándar simulada, 0.348 en R y 0.352 en Python, coincide con el valor de la fórmula $\sigma/\sqrt{S_{xx}} = 0.347$. Los dos lenguajes difieren solo en el último dígito porque usan generadores de números aleatorios distintos; la coincidencia con la teoría es lo importante.

El haz de rectas ajustadas de la Figure 15 se cree mucho mejor cuando lo construyes una muestra a la vez.

Qué observar. Toma muestras de una en una y ninguna $b_1$ es igual a 3.5702, y aun así el montón que construyen queda centrado ahí. Prueba esto. Presiona Tomar 1000 y compara el valor DE de $b_1$ con $\sigma/\sqrt{S_{xx}} = 0.3469$ que está a su lado: la fórmula de la varianza del Teorema 2.9 no es una aproximación, es el ancho de este montón. Volver a la Sección 2.5.

2.6 Máxima verosimilitud bajo errores normales¶

Intuición¶

Hasta ahora elegimos los mínimos cuadrados porque el error al cuadrado es una penalización razonable, y Gauss-Markov premió esa elección con la propiedad de mejor-lineal-insesgado, todo sin suponer ninguna forma para la distribución del error. Ahora supón que añadimos un supuesto: los errores están distribuidos normalmente. Un principio distinto, la máxima verosimilitud, pregunta qué valores de los parámetros hacen más probables los datos que realmente observamos. De forma notable, bajo errores normales ese principio apunta a la mismísima recta. La Figure 19 muestra los dos criterios, la suma de errores al cuadrado y el logaritmo de verosimilitud negativo, tocando fondo exactamente en la misma pendiente.

Fórmula¶

El modelo de regresión con errores normales (normal error regression model) (Definición 2.11) conserva todo lo de 2.1 El modelo de regresión lineal simple y añade una forma:

Bajo normalidad los errores no correlacionados ahora son plenamente independientes. La verosimilitud (likelihood) es la densidad conjunta de los $Y_i$ observados leída como función de los parámetros:

• $L$ es grande para valores de parámetros que hacen que los datos observados se vean típicos, pequeña para valores que los hacen ver raros.

• Las estimaciones de máxima verosimilitud son los $\beta_0, \beta_1, \sigma^2$ que maximizan $L$.

Derivación (la máxima verosimilitud es igual a los mínimos cuadrados)¶

Demostración. Maximizar $L$ es más fácil a través de su logaritmo, que se maximiza en el mismo lugar porque $\log$ es creciente:

Coeficientes. Fija cualquier $\sigma^2 > 0$. Los dos primeros términos no involucran a $\beta_0$ ni a $\beta_1$, y el último término tiene $-1/(2\sigma^2) < 0$ multiplicando a $\sum (Y_i - \beta_0 - \beta_1 X_i)^2$. Maximizar $\ln L$ sobre $\beta_0, \beta_1$ por lo tanto significa minimizar esa suma de cuadrados, que es exactamente el criterio de mínimos cuadrados $Q$ de 2.2 Mínimos cuadrados desde los primeros principios. Así que las estimaciones de máxima verosimilitud de los coeficientes son las estimaciones de mínimos cuadrados:

Varianza. Sustituye $b_0, b_1$ para que la suma de cuadrados se convierta en $\mathrm{SSE}$, luego deriva $\ln L$ respecto a $\sigma^2$ e iguala a cero:

La estimación de varianza de máxima verosimilitud divide entre $n$, no entre $n-2$. De 2.4 Estimar la varianza del error sabemos que el divisor $n$ es sesgado bajo, subestimando $\sigma^2$, porque ignora los dos grados de libertad gastados en $b_0$ y $b_1$. Por eso el libro, y toda la inferencia del Capítulo 3, usa el $\mathrm{MSE} = \mathrm{SSE}/(n-2)$ insesgado en lugar del $\hat\sigma^2_{\text{MLE}}$ de máxima verosimilitud. $\blacksquare$

R y Python¶

2.7 Resumen del capítulo¶

Este capítulo construyó una regresión lineal simple desde cero y demostró las propiedades que la hacen confiable. Partiendo del modelo $Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$ y sus tres supuestos sobre los errores, definimos el criterio de mínimos cuadrados y derivamos los estimadores de dos formas independientes, por cálculo a través de las ecuaciones normales y por una identidad algebraica que muestra que ninguna otra recta tiene una SSE menor. Los residuos obedecen identidades exactas que sirven a la vez como verificación de corrección; la varianza del error se estima con $\mathrm{MSE} = \mathrm{SSE}/(n-2)$, cuyo divisor $n-2$ es exactamente lo que la hace insesgada. Tratar $b_0$ y $b_1$ como aleatorios a lo largo de muestras repetidas dio sus medias y varianzas, el teorema de Gauss-Markov mostró que son los mejores entre los estimadores lineales insesgados sin necesidad de normalidad, y la máxima verosimilitud bajo errores normales reprodujo la mismísima recta.

Resultados clave de un vistazo

Para los datos de Toluca estos dan $b_0 = 62.37$, $b_1 = 3.5702$, $\mathrm{SSE} = 54{,}825$, $\mathrm{MSE} = 2384$, $s = 48.82$, $s\{b_1\} = 0.347$, $s\{b_0\} = 26.18$ y $\hat\sigma^2_{\text{MLE}} = 2193$.

Términos clave

Modelo de regresión lineal simple, función de regresión, respuesta media, error aleatorio, criterio de mínimos cuadrados, estimaciones de mínimos cuadrados, ecuaciones normales, valor ajustado, residuo, suma de cuadrados del error (SSE), cuadrado medio del error (MSE), homocedasticidad, estimador insesgado, teorema de Gauss-Markov, mejor estimador lineal insesgado (BLUE), modelo de regresión con errores normales, máxima verosimilitud.

Ahora deberías poder

• Enunciar el modelo de regresión lineal simple y explicar la función que cumple cada supuesto sobre el error.

• Derivar los estimadores de mínimos cuadrados de dos formas, por cálculo y por la identidad algebraica.

• Calcular $b_0$, $b_1$, los valores ajustados y los residuos a mano y con software.

• Demostrar las identidades de los residuos y usarlas como una verificación incorporada de un ajuste.

• Explicar el divisor $n-2$ derivando $E\{\mathrm{MSE}\} = \sigma^2$, y calcular $\mathrm{MSE}$ y $s$.

• Derivar las medias y varianzas de $b_0$ y $b_1$ y demostrar que $b_1$ es el mejor estimador lineal insesgado.

• Mostrar que la máxima verosimilitud bajo errores normales reproduce las estimaciones de mínimos cuadrados.

Dónde encaja esto. Este capítulo es la primera vuelta completa alrededor del flujo de trabajo de El flujo de trabajo del modelado. Llevó la pregunta de Toluca por PREGUNTAR y EXPLORAR, se asentó en AJUSTAR, y dio los primeros pasos hacia VERIFICAR (las identidades de los residuos y el gráfico de residuos) y USAR (leer la pendiente en horas por unidad). El Capítulo 1 argumentó que la regresión es una forma de pensar sobre cómo una cantidad se mueve con otra; este capítulo convirtió esa idea en una recta específica y calculable y demostró las propiedades que la hacen confiable. Todo aquí es estimación puntual: conjeturas únicas y mejores para $\beta_0$, $\beta_1$ y $\sigma^2$. Lo que todavía no hemos hecho es cuantificar la incertidumbre con intervalos y pruebas. Ese es el Capítulo 3, que toma las distribuciones muestrales bosquejadas aquí, añade el supuesto de errores normales por completo, y construye intervalos $t$ para la pendiente y el intercepto, la descomposición ANOVA (3.6 El enfoque del análisis de varianza), la prueba $F$ (3.7 La prueba F y su equivalencia con la prueba t), y la diferencia entre un intervalo de confianza para una respuesta media y un intervalo de predicción para una corrida nueva (3.5 Intervalo de predicción para una observación nueva). El ajuste exacto de Toluca fijado aquí ($b_1 = 3.5702$, $\mathrm{MSE} = 2384$) se traslada sin cambios, ahora con inferencia adjunta. Viaja aún más lejos: el Capítulo 5 vuelve a pasar la misma pendiente por una prueba de permutación (5.2 La prueba de permutación para la pendiente) y compara su valor p con el del Capítulo 3, y en el Capítulo 7 los mismos estimadores reaparecen en forma matricial (7.1 El modelo y los mínimos cuadrados en forma matricial), donde las ecuaciones normales de dos líneas se convierten en una sola ecuación $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$ y la demostración de Gauss-Markov (7.6 El teorema de Gauss-Markov) se generaliza a cualquier número de predictores a la vez.

2.8 Preguntas frecuentes¶

P1. ¿Por qué elevar al cuadrado los residuos en lugar de tomar valores absolutos? Ambas son pérdidas razonables. Elevar al cuadrado hace que el objetivo sea suave y diferenciable, así que el cálculo da una solución limpia en forma cerrada (las ecuaciones normales), y es exactamente la pérdida para la cual se cumplen los resultados de Gauss-Markov y de verosimilitud normal. La pérdida de valor absoluto da la regresión de “mínimas desviaciones absolutas”, que es más resistente a los valores atípicos pero no tiene fórmula simple y necesita resolverse iterativamente. El curso empieza con el error al cuadrado porque es donde la teoría es más limpia.

P2. ¿El intercepto $b_0 = 62.37$ es realmente “horas para producir cero unidades”? No, no de forma significativa. Un tamaño de lote de cero está muy fuera de los datos (la corrida más pequeña fue de 20 unidades), así que el intercepto es una extrapolación. Léelo como el número que ancla la altura de la recta, calculado para que la recta pase por $(\bar{X}, \bar{Y})$, no como un tiempo de preparación físico. Interpretar un intercepto literalmente solo tiene sentido cuando $X = 0$ está dentro del rango observado o cerca de él.

P3. ¿Cuál es la diferencia entre un error $\varepsilon_i$ y un residuo $e_i$? El error $\varepsilon_i = Y_i - (\beta_0 + \beta_1 X_i)$ es la brecha a la recta verdadera y desconocida, y nunca podemos observarlo. El residuo $e_i = Y_i - \hat{Y}_i$ es la brecha a nuestra recta estimada, y lo calculamos directamente. Los residuos son nuestras estimaciones visibles de los errores invisibles, y por eso estimamos $\sigma^2$ a partir de ellos.

P4. ¿Por qué dividir la SSE entre $n-2$ y no entre $n-1$ como una varianza muestral ordinaria? Una varianza muestral divide entre $n-1$ porque estima una cantidad, la media, a partir de los datos. La regresión estima dos, el intercepto y la pendiente, antes de medir los residuos, así que gasta dos grados de libertad, dejando $n-2$. La derivación en 2.4 Estimar la varianza del error muestra $E\{\mathrm{SSE}\} = (n-2)\sigma^2$ exactamente, así que $n-2$ es el divisor que hace insesgado a $\mathrm{MSE}$.

P5. ¿Necesito que los errores sean normales para que los mínimos cuadrados funcionen? No. Las estimaciones de mínimos cuadrados, su insesgadez, sus fórmulas de varianza y la optimalidad de Gauss-Markov se cumplen todas sin ningún supuesto sobre la forma de la distribución del error, solo media cero, varianza constante y errores no correlacionados. La normalidad es un supuesto extra, añadido en 2.6 Máxima verosimilitud bajo errores normales y usado en el Capítulo 3 para construir intervalos y pruebas $t$ y $F$ exactos.

P6. ¿Por qué el haz de rectas ajustadas en la Figure 15 se estrecha en el medio? Toda recta de mínimos cuadrados pasa por $(\bar{X}, \bar{Y})$ (propiedad 5 de los residuos), así que cerca de $\bar{X}$ las rectas solo pueden temblar un poco. Lejos de $\bar{X}$, un pequeño cambio en la pendiente hace pivotar la recta una gran distancia vertical, así que las rectas se abren. La varianza de un valor ajustado crece con la distancia desde $\bar{X}$, lo cual el Capítulo 3 convierte en la forma curva de una banda de confianza.

P7. R imprime 4.45e-10 y Python imprime 4.45e-10 para el valor p de la pendiente, pero no supusimos normalidad. ¿De dónde vino eso? Esos valores p y los valores t de los errores estándar sí usan el modelo de errores normales, y los justificaremos apropiadamente en el Capítulo 3. Este capítulo muestra los errores estándar en sí, que vienen de las fórmulas de varianza y no necesitan normalidad; las estrellas de significancia y los valores p son un adelanto de la inferencia que viene.

2.9 Problemas de práctica¶

1. (A) Enuncia el modelo de regresión lineal simple y lista sus tres supuestos sobre los errores, diciendo en una frase qué garantiza cada supuesto.

2. (A) Explica la diferencia entre la función de regresión $E\{Y\} = \beta_0 + \beta_1 X$ y una recta ajustada $\hat{Y} = b_0 + b_1 X$. ¿Cuál involucra parámetros y cuál involucra estimaciones?

3. (A) En el ajuste de Toluca, interpreta la pendiente 3.5702 en las unidades del problema. ¿Por qué es inseguro usar la recta para predecir horas para un lote de 500 unidades?

4. (A) Da una interpretación de una oración de $s = 48.82$ horas para un gerente que nunca ha tomado estadística.

5. (A) ¿Por qué el intercepto de mínimos cuadrados para los datos de Toluca no es un “tiempo de preparación” significativo? ¿Qué tendría que ser cierto de los datos para que un intercepto fuera interpretable?

6. (A) Sin calcular, enuncia los valores de $\sum e_i$ y $\sum X_i e_i$ para cualquier ajuste de mínimos cuadrados con intercepto, y di qué ecuación normal da cada uno.

7. (A) Una recta ajustada tiene $\sum e_i = 0$ pero $\sum X_i e_i = 120$. ¿Puede ser un ajuste de mínimos cuadrados? Explica.

8. (A) Explica por qué dispersar los valores del predictor más lejos entre sí (aumentar $S_{xx}$) hace más precisa la estimación de la pendiente, refiriéndote a su fórmula de varianza.

9. (B) Partiendo de $Q(b_0, b_1) = \sum (Y_i - b_0 - b_1 X_i)^2$, deriva ambas ecuaciones normales por derivación, y resuélvelas para $b_0$ y $b_1$ (Teorema 2.4).

10. (B) Demuestra la identidad algebraica $\sum (Y_i - a - cX_i)^2 = \mathrm{SSE} + \sum g_i^2$ para una recta arbitraria $(a, c)$, identificando $g_i$ y mostrando que el término cruzado es cero. Concluye que los mínimos cuadrados son el minimizador único.

11. (B) Muestra que la recta de mínimos cuadrados pasa por $(\bar{X}, \bar{Y})$, y que $\sum \hat{Y}_i = \sum Y_i$ (Teorema 2.6).

12. (B) Escribiendo $b_1 = \sum k_i Y_i$ con $k_i = (X_i - \bar{X})/S_{xx}$, demuestra $\sum k_i = 0$, $\sum k_i X_i = 1$ y $\sum k_i^2 = 1/S_{xx}$, luego úsalos para mostrar $E\{b_1\} = \beta_1$ y $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$.

13. (B) Deriva $\operatorname{Var}\{b_0\} = \sigma^2(1/n + \bar{X}^2/S_{xx})$ escribiendo $b_0 = \sum m_i Y_i$ con $m_i = 1/n - \bar{X}k_i$.

14. (B) Demuestra el teorema de Gauss-Markov para la pendiente (Teorema 2.10): cualquier estimador lineal insesgado $\sum c_i Y_i$ de $\beta_1$ tiene varianza al menos $\sigma^2/S_{xx}$, con igualdad solo para $b_1$.

15. (B) Deriva $E\{\mathrm{SSE}\} = (n-2)\sigma^2$ (Teorema 2.8) usando $\mathrm{SSE} = S_{yy} - b_1^2 S_{xx}$ y los hechos $E\{b_1\} = \beta_1$, $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$.

16. (B) Escribe el logaritmo de verosimilitud para el modelo con errores normales y muestra que maximizarlo sobre $\beta_0, \beta_1$ es equivalente a los mínimos cuadrados (Teorema 2.12). Luego deriva $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$.

17. (B) Muestra que $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$ es sesgado para $\sigma^2$, y calcula el sesgo exacto en términos de $\sigma^2$ y $n$.

18. (B) Demuestra la descomposición $\mathrm{SSE} = S_{yy} - b_1^2 S_{xx}$ partiendo de $e_i = (Y_i - \bar{Y}) - b_1(X_i - \bar{X})$.

19. (B) Un estudiante propone estimar la pendiente por el cociente de medias $\bar{Y}/\bar{X}$ (una recta por el origen). Muestra que esto es en general sesgado para $\beta_1$, y explica qué supuesto del modelo impone implícita, y erróneamente.

20. (C) Lee toluca.csv en R o Python, ajusta el modelo, y reproduce $b_0$, $b_1$ y $s$. Confirma $b_1 = S_{xy}/S_{xx}$ a partir de las sumas crudas.

21. (C) Calcula los residuos en software y verifica numéricamente que $\sum e_i = 0$, $\sum X_i e_i = 0$ y $\sum \hat{Y}_i = \sum Y_i$. Reporta las tres sumas.

22. (C) Usando toluca_mini.csv, calcula $b_0$ y $b_1$ a mano (muestra la tabla de desviaciones), luego reprodúcelos en software. Explica por qué difieren de las estimaciones de datos completos.

23. (C) Predice las horas de trabajo medias para tamaños de lote 40, 80 y 120 usando la recta ajustada, y marca cuáles de estos son interpolación y cuáles (si alguno) son extrapolación.

24. (C) Ajusta el modelo, extrae $\mathrm{MSE}$ de dos formas (de summary / fit.summary() y a mano como $\mathrm{SSE}/(n-2)$), y confirma que coinciden. Luego calcula $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$ y reporta la diferencia porcentual.

25. (C) Escribe una simulación corta (semilla 4210) que extraiga 2000 muestras del modelo ajustado en los tamaños de lote de Toluca, reajuste cada una, y reporte la media y la desviación estándar de las 2000 estimaciones de intercepto $b_0$. Compara la desviación estándar con la fórmula $s\{b_0\}$.

26. (C) Haz el gráfico de residuos contra tamaño de lote en R o Python. Describe cómo se vería una violación del supuesto de varianza constante en este gráfico, y di si los residuos de Toluca la muestran.

27. (C) Invierte los papeles: regresa lotsize sobre hours. Muestra que la nueva pendiente no es el recíproco de la pendiente original, y explica por qué regresar $Y$ sobre $X$ y $X$ sobre $Y$ dan rectas distintas.

28. (B) Un analista tiene solo los números de resumen $\bar{X} = 70$, $\bar{Y} = 312.28$, $S_{xx} = 19{,}800$ y $S_{xy} = 70{,}690$, con los datos crudos perdidos. ¿Pueden recuperar $b_0$ y $b_1$? ¿Pueden recuperar $\mathrm{MSE}$? Explica qué es y qué no es reconstruible a partir de estos cuatro números.