3. Inferencia para la regresión lineal simple (en espanol)

3. Inferencia para la regresión lineal simple¶

En 2.2 Mínimos cuadrados desde los primeros principios ajustamos una recta a las 25 corridas de producción de la Toluca Company y obtuvimos una pendiente de $b_1 = 3.5702$ horas por unidad de tamaño de lote. Los ingenieros quedaron satisfechos y enseguida hicieron las preguntas que cualquier persona cuidadosa hace sobre un número: ¿qué tan seguros están? Otras 25 corridas habrían dado una pendiente algo distinta. ¿Podría la pendiente verdadera ser realmente cero, es decir, que el tamaño del lote no afecte en nada a las horas y nos estemos engañando? Y la pregunta práctica, la razón por la que recogieron los datos: si programamos una nueva corrida de 100 unidades la próxima semana, ¿cuántas horas debemos presupuestar y qué tan lejos podríamos quedar?

El Capítulo 2 dio estimaciones puntuales: una sola mejor conjetura para el intercepto, la pendiente y la varianza del error. Una estimación puntual sin ninguna medida de su incertidumbre es media respuesta. Este capítulo aporta la otra mitad. Adjuntamos márgenes de error a la pendiente y al intercepto, probamos si la pendiente es realmente distinta de cero, y, lo más importante para la planeación, construimos dos intervalos muy diferentes alrededor de una predicción. Uno dice dónde cae la corrida promedio de un tamaño dado. El otro dice dónde cae una sola corrida nueva. La Figure 1 muestra ambos trazados alrededor de la recta de Toluca, y la brecha entre ellos es la idea más útil del capítulo.

Todo lo que sigue descansa sobre un supuesto adicional que el Capítulo 2 introdujo pero aún no había usado para la inferencia: los errores están distribuidos normalmente (2.6 Máxima verosimilitud bajo errores normales). Con ese supuesto las distribuciones muestrales bosquejadas en 2.5 Comportamiento muestral y el teorema de Gauss-Markov se vuelven distribuciones $t$ y $F$ exactas, y de ahí se siguen intervalos y pruebas honestos.

3.1 Las distribuciones muestrales de $b_0$ y $b_1$¶

Intuición¶

Ya sabemos por 2.5 Comportamiento muestral y el teorema de Gauss-Markov que la estimación de la pendiente $b_1$ es insesgada, con varianza $\sigma^2/S_{xx}$, donde $S_{xx} = \sum_i (X_i - \bar{X})^2$ es la dispersión de los valores del predictor alrededor de su media $\bar{X}$. Eso nos dice dónde se centra $b_1$ y cuánto se dispersa, pero no la forma de su distribución. Para construir un intervalo necesitamos la forma. El truco es que $b_1$ es una suma ponderada de las respuestas, y una suma ponderada de variables aleatorias normales independientes es también normal. Así que, en cuanto suponemos errores normales, $b_1$ es exactamente normal, y podemos estandarizarla igual que la estadística introductoria estandariza una media muestral.

Hay una complicación. Estandarizar necesita la desviación estándar verdadera del error $\sigma$, que no conocemos. La reemplazamos por la estimación $s = \sqrt{\mathrm{MSE}}$. Cambiar una constante conocida por una estimación ruidosa añade un poco de incertidumbre extra, y esa incertidumbre extra es exactamente lo que convierte la distribución normal en una distribución $t$. La Figure 2 lo confirma: estandariza 5000 pendientes simuladas por sus propios errores estándar estimados, y el montón coincide con una curva $t$, con colas un poco más pesadas que la normal.

Fórmula¶

Bajo el modelo de errores normales $Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$ con $\varepsilon_i \overset{\text{iid}}{\sim} N(0, \sigma^2)$, los estimadores son exactamente normales:

En palabras: cada estimador está centrado en el parámetro que estima, con la varianza que derivamos en el Capítulo 2, y ahora con una forma normal. Para estandarizar cualquiera de los dos necesitamos la dispersión de la estimación. El error estándar estimado de la pendiente (Definición 3.1) la proporciona.

El error estándar estimado $s\{b_1\}$ es nuestra mejor conjetura de cuánto rebota $b_1$ de una muestra a la siguiente. Estandarizar $b_1$ con él da un estadístico $t$, el resultado que se enuncia a continuación.

En palabras: resta el valor verdadero y divide por el error estándar estimado, y el resultado sigue una distribución $t$ con $n-2$ grados de libertad. Los grados de libertad son $n-2$ porque ese es el divisor en $\mathrm{MSE}$, y los dos grados perdidos pagan por $b_0$ y $b_1$ (2.4 Estimar la varianza del error). Este único resultado $t$ alimenta cada intervalo y prueba en el resto del capítulo.

Esa última frase merece una imagen. Casi todo lo que viene es una sola idea, el estadístico $t$ de arriba, apuntado a una estimación distinta cada vez. La Figure 3 es el mapa del capítulo en una sola página: tenla presente y cada sección nueva se vuelve “el mismo movimiento, nuevo objetivo”.

Derivación¶

Derivación (el estadístico $t$ para la pendiente). Armamos el resultado a partir de tres hechos.

Normalidad. De 2.5 Comportamiento muestral y el teorema de Gauss-Markov, $b_1 = \sum_i k_i Y_i$ con $k_i = (X_i - \bar{X})/S_{xx}$. Bajo el modelo de errores normales los $Y_i$ son variables normales independientes, y cualquier combinación lineal de normales independientes es normal. Con la media y la varianza ya derivadas, $b_1 \sim N(\beta_1, \sigma^2/S_{xx})$. Dividir por la desviación estándar verdadera da una normal estándar:

Una chi cuadrada, independiente de la pendiente. Se puede demostrar que

y que esta cantidad es estadísticamente independiente de $b_1$ (de hecho de todo el vector $(b_0, b_1)$). La demostración limpia usa la matriz sombrero y un teorema sobre formas cuadráticas en vectores normales; la damos completa en 7.3 La matriz sombrero, una vez que tengamos las herramientas matriciales. Por ahora tomamos estos dos hechos por fe, y notamos que la simulación en la Figure 2 es una comprobación directa de que el resultado final es correcto.

Combinar en una $t$. Una normal estándar dividida por la raíz cuadrada de una chi cuadrada independiente sobre sus grados de libertad es, por definición, una variable aleatoria $t$:

El paso intermedio solo desarrolla el numerador y el denominador: el numerador es $Z$, y $s\{b_1\}/(\sigma/\sqrt{S_{xx}}) = \sqrt{\mathrm{MSE}/\sigma^2} = \sqrt{\chi^2_{n-2}/(n-2)}$. El argumento idéntico con los pesos $m_i$ de 2.5 Comportamiento muestral y el teorema de Gauss-Markov da el estadístico $t$ para $b_0$. $\blacksquare$

El código de preparación de abajo reajusta el modelo de Toluca y recupera los errores estándar que usaremos en todo el capítulo.

toluca <- read.csv("data/toluca.csv")
fit <- lm(hours ~ lotsize, data = toluca)
n <- nrow(toluca)
x <- toluca$lotsize
y <- toluca$hours
Sxx <- sum((x - mean(x))^2)
MSE <- sum(residuals(fit)^2) / (n - 2)
round(c(n = n, Sxx = Sxx, MSE = MSE, s = sqrt(MSE),
        se_b1 = sqrt(MSE / Sxx)), 4)

         n        Sxx        MSE          s      se_b1 
   25.0000 19800.0000  2383.7156    48.8233     0.3470 

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

toluca = pd.read_csv("data/toluca.csv")
fit = smf.ols("hours ~ lotsize", data=toluca).fit()
n = len(toluca)
x = toluca["lotsize"].to_numpy()
y = toluca["hours"].to_numpy()
Sxx = np.sum((x - x.mean()) ** 2)
MSE = np.sum(fit.resid ** 2) / (n - 2)
print(round(MSE, 4), round(np.sqrt(MSE), 4), round(np.sqrt(MSE / Sxx), 4))

2383.7156 48.8233 0.347

El error estándar estimado de la pendiente es $s\{b_1\} = 0.347$ horas por unidad, el mismo número que el Capítulo 2 calculó a partir de la fórmula de la varianza. Es la vara de medir para todo lo que sigue.

3.2 Inferencia para la pendiente¶

Intuición¶

La pendiente suele ser el número del que trata un estudio: ¿el tamaño del lote afecta a las horas, y por cuánto? Dos preguntas vienen en pareja. Un intervalo de confianza (Definición 3.3) responde “¿por cuánto?” al dar un rango de valores plausibles para $\beta_1$. Una prueba de hipótesis responde “¿hay algún efecto en absoluto?” al preguntar si cero es un valor plausible. Ambas se construyen a partir del único estadístico $t$ de 3.1 Las distribuciones muestrales de $b_0$ y $b_1$, así que siempre concuerdan: un intervalo del 95% excluye a cero exactamente cuando la prueba bilateral rechaza al nivel del 5%.

Fórmula¶

En palabras: toma la estimación y extiéndete un multiplicador de errores estándar a cada lado. Para probar $H_0: \beta_1 = 0$ contra $H_a: \beta_1 \ne 0$, forma el estadístico de prueba

y rechaza $H_0$ cuando $|t^{\ast}| > t_{1-\alpha/2,\,n-2}$, o equivalentemente cuando el valor p bilateral $2\,P(t_{n-2} > |t^{\ast}|)$ está por debajo de $\alpha$. En palabras: si la estimación queda a muchos errores estándar de cero, cero no es creíble. Nada aquí es especial de cero; para probar $H_0: \beta_1 = c$ para cualquier valor $c$, pon $c$ en lugar del 0.

El 95% de “95% de confianza” es una promesa sobre muestras repetidas, no sobre el único intervalo que tienes enfrente, y eso cuesta aceptarlo de palabra. El simulador de abajo toma las muestras para que veas cómo se cumple la promesa y, de vez en cuando, cómo se rompe.

3.3 Inferencia para el intercepto¶

El intercepto es donde la recta ajustada cruza el eje vertical: las horas que el modelo predice para un lote de cero unidades. Es el mismo tipo de estimación que la pendiente, así que recibe exactamente la misma maquinaria. Sigue el mapa de Figure 3 una casilla más allá, y solo cambia el error estándar.

El intercepto recibe el mismo tratamiento, con su propio error estándar $s\{b_0\} = \sqrt{\mathrm{MSE}\,(1/n + \bar{X}^2/S_{xx})}$. El intervalo de confianza es $b_0 \pm t_{1-\alpha/2,\,n-2}\, s\{b_0\}$, y la prueba de $H_0: \beta_0 = 0$ usa $t^{\ast} = b_0/s\{b_0\}$.

La razón por la que el intercepto va en segundo lugar es que suele ser menos interesante y a menudo no interpretable. Para Toluca, $\beta_0$ es las horas medias en un tamaño de lote de cero, una extrapolación muy fuera de los datos (2.2 Mínimos cuadrados desde los primeros principios), así que una prueba de $\beta_0 = 0$ responde una pregunta que nadie hizo. El intervalo aún vale la pena calcularlo, aunque solo sea para ver qué tan flojamente queda sujeto el intercepto.

confint(fit, "(Intercept)", level = 0.95)

               2.5 %  97.5 %
(Intercept) 8.213711 116.518

print(fit.conf_int(alpha=0.05).loc["Intercept"])

0      8.213711
1    116.518006
Name: Intercept, dtype: float64

El intervalo del 95% para el intercepto va de unas 8 a 117 horas, un tramo de más de 100 horas. Ese ancho enorme es el precio de extrapolar a $X = 0$, que queda 70 unidades por debajo del lote más pequeño observado. Compáralo con el intervalo ajustado de la pendiente: la pendiente está anclada por datos a ambos lados, el intercepto no.

3.4 Intervalo de confianza para una respuesta media¶

Intuición¶

Ahora pasamos de los coeficientes a una predicción. Pregunta: entre todas las corridas de tamaño de lote $X_h = 65$, ¿cuál es el número promedio de horas? Ese promedio es el punto sobre la recta verdadera, $E\{Y_h\} = \beta_0 + \beta_1 X_h$, y nuestra estimación de él es el valor ajustado $\hat{Y}_h = b_0 + b_1 X_h$. Como $b_0$ y $b_1$ oscilan de una muestra a otra, también lo hace $\hat{Y}_h$, y queremos un intervalo para la media verdadera.

La característica clave es que el ancho del intervalo depende de dónde se ubica $X_h$. Cerca del centro $\bar{X}$ la recta ajustada queda sujeta con firmeza, así que el intervalo es angosto. Lejos de $\bar{X}$ una pequeña inclinación en la pendiente hace subir o bajar mucho la recta, así que el intervalo se ensancha. Este es el estrechamiento que vimos en el haz de rectas ajustadas de Figure 15, ahora convertido en una fórmula.

Fórmula¶

La varianza estimada de la media ajustada en $X_h$ es

y de ahí se sigue el intervalo de confianza para la respuesta media.

En palabras: cuanto más lejos está $X_h$ del valor promedio del predictor $\bar{X}$, mayor es el segundo término dentro del paréntesis, y más ancho el intervalo. En $X_h = \bar{X}$ el segundo término se anula y el intervalo está en su punto más angosto, con el ancho fijado solo por $1/n$.

Derivación¶

Derivación (varianza de la media ajustada). Escribe el valor ajustado en su forma pivoteada, usando $b_0 = \bar{Y} - b_1 \bar{X}$:

Las dos piezas $\bar{Y}$ y $b_1$ no están correlacionadas. Para verlo, recuerda $b_1 = \sum k_i Y_i$ y $\bar{Y} = \sum \tfrac{1}{n} Y_i$, así que, usando la independencia de los $Y_i$ y $\sum k_i = 0$,

Como las dos piezas no están correlacionadas, sus varianzas se suman:

Reemplazar $\sigma^2$ por $\mathrm{MSE}$ da $s^2\{\hat{Y}_h\}$. El intervalo $t$ se sigue entonces del mismo argumento de estandarizar y dividir de 3.1 Las distribuciones muestrales de $b_0$ y $b_1$, ya que $(\hat{Y}_h - E\{Y_h\})/s\{\hat{Y}_h\} \sim t_{n-2}$. $\blacksquare$

La Figure 5 grafica $s\{\hat{Y}_h\}$ contra $X_h$: una curva que baja a su mínimo en $\bar{X} = 70$ y sube simétricamente a ambos lados. Esa curva es lo que le da a la banda de confianza de Figure 1 su suave forma de moño.

3.5 Intervalo de predicción para una observación nueva¶

Intuición¶

Aquí está la distinción que confunde a más profesionales que ninguna otra. El intervalo de confianza de 3.4 Intervalo de confianza para una respuesta media trata sobre la corrida promedio de tamaño $X_h$. Pero los ingenieros no programan una corrida promedio; programan una corrida específica la próxima semana, y quieren saber dónde caerán sus horas. Una sola corrida se dispersa alrededor de la media por su propio error aleatorio $\varepsilon$, y ninguna cantidad de datos elimina esa dispersión. Así que predecir una sola observación nueva conlleva dos fuentes de incertidumbre: nuestra incertidumbre sobre dónde está la recta (el mismo $s^2\{\hat{Y}_h\}$ de antes), más la propia varianza irreducible de la corrida nueva $\sigma^2$. El intervalo que contempla ambas es el intervalo de predicción (Definición 3.6), y siempre es más ancho que el intervalo de confianza.

Toda la elección se reduce a una pregunta, y la Figure 6 la expone. Pregúntate qué estás acotando: el resultado promedio sobre muchos casos en $X_h$, o el resultado de un caso nuevo específico. La respuesta de la izquierda es un intervalo de confianza; la de la derecha es un intervalo de predicción. Decide la pregunta primero, y la fórmula se sigue.

Fórmula¶

Para una observación nueva $Y_{h(\text{new})}$ en $X_h$, la varianza estimada del error de predicción es $s^2\{\text{pred}\} = \mathrm{MSE}(1 + 1/n + (X_h - \bar{X})^2/S_{xx})$, y de ahí se sigue el intervalo de predicción.

En palabras: esta es la varianza de la respuesta media con un $\mathrm{MSE}$ extra (el 1 inicial) añadido por la propia dispersión de la corrida nueva. Ese término extra suele ser el dominante, por lo cual el intervalo de predicción es mucho más ancho y por lo cual, a diferencia del intervalo de confianza, nunca se encoge a cero ni siquiera con datos infinitos.

Derivación¶

Derivación (varianza de una predicción). La cantidad que acotamos es el error de predicción $Y_{h(\text{new})} - \hat{Y}_h$, la brecha entre la corrida nueva real y nuestro pronóstico. La corrida nueva es $Y_{h(\text{new})} = \beta_0 + \beta_1 X_h + \varepsilon_{\text{new}}$, y su error $\varepsilon_{\text{new}}$ es independiente de los datos pasados que produjeron $\hat{Y}_h$. Dos piezas independientes, así que sus varianzas se suman:

Reemplazar $\sigma^2$ por $\mathrm{MSE}$ da $s^2\{\text{pred}\}$, y estandarizar el error de predicción por $s\{\text{pred}\}$ produce de nuevo una distribución $t_{n-2}$. El intervalo $t$ se sigue. $\blacksquare$

La brecha entre los dos intervalos no es una elección de modelado; es un hecho sobre las dos preguntas, y cada intervalo mantiene su propia promesa del 95%. Podemos comprobarlo por simulación, que es la manera más segura de ver que un intervalo de predicción realmente atrapa una observación nueva el 95% de las veces. La Figure 7 corre el experimento.

set.seed(4210)
b0_true <- 62.366; b1_true <- 3.5702; sigma <- 48.82
xh <- 100
mu_h <- b0_true + b1_true * xh
tcrit <- qt(0.975, n - 2)
ci_hits <- 0; pi_hits <- 0; N <- 5000
for (i in 1:N) {
  ysim <- b0_true + b1_true * x + rnorm(n, 0, sigma)
  f <- lm(ysim ~ x)
  pr <- predict(f, data.frame(x = xh), se.fit = TRUE)
  s_mean <- pr$se.fit
  s_pred <- sqrt(pr$se.fit^2 + sum(residuals(f)^2) / (n - 2))
  if (abs(pr$fit - mu_h) <= tcrit * s_mean) ci_hits <- ci_hits + 1
  ynew <- mu_h + rnorm(1, 0, sigma)
  if (abs(pr$fit - ynew) <= tcrit * s_pred) pi_hits <- pi_hits + 1
}
round(c(CI_covers_mean = ci_hits / N, PI_covers_new = pi_hits / N), 4)

CI_covers_mean  PI_covers_new 
        0.9496         0.9458 

from scipy import stats

rng = np.random.default_rng(4210)
b0_true, b1_true, sigma = 62.366, 3.5702, 48.82
xh = 100.0
mu_h = b0_true + b1_true * xh
tcrit = stats.t.ppf(0.975, n - 2)
ci_hits = pi_hits = 0
N = 5000
for _ in range(N):
    ysim = b0_true + b1_true * x + rng.normal(0, sigma, size=n)
    b1s = np.sum((x - x.mean()) * (ysim - ysim.mean())) / Sxx
    b0s = ysim.mean() - b1s * x.mean()
    yhat = b0s + b1s * xh
    mse = np.sum((ysim - (b0s + b1s * x)) ** 2) / (n - 2)
    s_mean = np.sqrt(mse * (1 / n + (xh - x.mean()) ** 2 / Sxx))
    s_pred = np.sqrt(mse * (1 + 1 / n + (xh - x.mean()) ** 2 / Sxx))
    if abs(yhat - mu_h) <= tcrit * s_mean:
        ci_hits += 1
    ynew = mu_h + rng.normal(0, sigma)
    if abs(yhat - ynew) <= tcrit * s_pred:
        pi_hits += 1
print(round(ci_hits / N, 4), round(pi_hits / N, 4))

0.9492 0.9544

A lo largo de 5000 conjuntos de datos simulados, el intervalo de confianza cubre la media verdadera cerca del 95% de las veces y el intervalo de predicción cubre una corrida nueva cerca del 95% de las veces. Cada intervalo es honesto sobre su propia pregunta. El intervalo de predicción se gana su ancho extra; no está siendo cauteloso por gusto.

La separación entre los dos intervalos se cree mejor cuando puedes abrirla y cerrarla tú mismo. El simulador de abajo pone el tamaño de muestra y el punto de predicción bajo tu dedo, para que descubras cuál banda responde a más datos y cuál no.

3.6 El enfoque del análisis de varianza¶

Intuición¶

Hay una segunda manera, equivalente, de ver la misma inferencia, y escala muy bien hacia la regresión múltiple de capítulos posteriores. La idea es contabilidad de la variación. Los valores de la respuesta $Y_i$ varían; parte de esa variación la explica la recta, y parte la deja como dispersión residual. Si podemos partir la variación total limpiamente en una parte “explicada” y una parte “sobrante”, podemos compararlas, y una parte explicada grande respecto a la sobrante es evidencia de que la pendiente es real.

La Figure 9 muestra la partición para una sola corrida: la desviación total respecto a la media, $Y_i - \bar{Y}$, se descompone en la parte que la recta explica, $\hat{Y}_i - \bar{Y}$, más el residuo, $Y_i - \hat{Y}_i$. Notablemente, cuando elevas al cuadrado y sumas sobre todas las corridas, el término cruzado se anula y los totales se parten con la misma limpieza.

Fórmula¶

Define las tres sumas de cuadrados:

En palabras: la variación total de $Y$ respecto a su media (SSTO) es igual a la variación que la regresión explica (SSR, la dispersión de los valores ajustados respecto a la media) más la variación que no explica (SSE, la suma de cuadrados del error de 2.4 Estimar la varianza del error). Cada suma tiene grados de libertad: $n-1$ para SSTO, 1 para SSR (un predictor), y $n-2$ para SSE. Los grados de libertad se suman de la misma manera, $n-1 = 1 + (n-2)$. Dividir cada suma de cuadrados por sus grados de libertad da un cuadrado medio (Definición 3.8).

En palabras: un cuadrado medio es una suma de cuadrados promediada sobre sus grados de libertad, así que es un número tipo varianza. MSR y MSE son los dos que pondremos en un cociente para probar la pendiente.

Derivación¶

Derivación ($\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$). Parte cada desviación total encaminándola a través del valor ajustado:

Eleva ambos lados al cuadrado y suma sobre $i$:

El término cruzado es cero. Escribe $Y_i - \hat{Y}_i = e_i$ y desarrolla:

Ambas piezas se anulan por las identidades de residuos de 2.3 Valores ajustados y las propiedades de los residuos: $\sum \hat{Y}_i e_i = 0$ (los valores ajustados son ortogonales a los residuos) y $\sum e_i = 0$. Así que el término cruzado es cero y $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$. Las identidades de residuos que probamos en el Capítulo 2 vuelven a hacer el trabajo. $\blacksquare$

Derivación (cuadrados medios esperados). De 2.4 Estimar la varianza del error ya tenemos $E\{\mathrm{MSE}\} = \sigma^2$. Un cálculo de esperanza similar, usando $\mathrm{SSR} = b_1^2 S_{xx}$ y $E\{b_1^2\} = \operatorname{Var}\{b_1\} + \beta_1^2 = \sigma^2/S_{xx} + \beta_1^2$, da

En palabras: MSE siempre estima $\sigma^2$, mientras que MSR estima $\sigma^2$ más un término que es positivo exactamente cuando $\beta_1 \ne 0$. Así que si la pendiente es cero, MSR y MSE estiman lo mismo y su cociente debería estar cerca de 1; si la pendiente es distinta de cero, MSR está inflado y el cociente es grande. Ese cociente es la prueba $F$. $\blacksquare$

La tabla ANOVA reúne estas cantidades en un formato estándar. Este es el anclaje que capítulos posteriores extienden a sumas de cuadrados adicionales (8.3 Sumas de cuadrados adicionales).

3.7 La prueba F y su equivalencia con la prueba t¶

Intuición¶

La tabla ANOVA nos entrega una prueba casi gratis. Bajo $H_0: \beta_1 = 0$, tanto MSR como MSE estiman el mismo $\sigma^2$, así que su cociente $F^{\ast} = \mathrm{MSR}/\mathrm{MSE}$ ronda 1. Bajo $H_a: \beta_1 \ne 0$, MSR está inflado por $\beta_1^2 S_{xx}$ y el cociente crece. Un $F^{\ast}$ grande es evidencia contra la nula.

Pero ya probamos $\beta_1 = 0$ con una prueba $t$ en 3.2 Inferencia para la pendiente y obtuvimos $t^{\ast} = 10.29$. ¿Concuerdan las dos pruebas? Sí, exacta y siempre: en la regresión lineal simple el estadístico $F$ es el cuadrado del estadístico $t$. Esto vale la pena probarlo, porque muestra que la prueba $F$ de ANOVA no es una prueba nueva, solo la prueba $t$ con otra ropa.

Fórmula¶

La prueba $F$ de $H_0: \beta_1 = 0$ usa

y rechaza $H_0$ cuando $F^{\ast} > F_{1-\alpha;\,1,\,n-2}$. La identidad afirmada es

donde $t^{\ast} = b_1/s\{b_1\}$ es el estadístico $t$ de la pendiente. En palabras: el estadístico $F$ y el estadístico $t$ al cuadrado son el mismo número, y las pruebas dan valores p idénticos.

Derivación¶

Derivación ($F^{\ast} = (t^{\ast})^2$). Parte de la suma de cuadrados de la regresión. Como $\hat{Y}_i - \bar{Y} = b_1 (X_i - \bar{X})$,

Con un grado de libertad de regresión, $\mathrm{MSR} = \mathrm{SSR} = b_1^2 S_{xx}$. Por lo tanto

Ahora escribe el estadístico $t$ y elévalo al cuadrado, usando $s^2\{b_1\} = \mathrm{MSE}/S_{xx}$:

Las dos expresiones son literalmente la misma. Distribucionalmente esto coincide con un hecho estándar: el cuadrado de una variable aleatoria $t_{n-2}$ es una variable aleatoria $F_{1,\,n-2}$, así que elevar al cuadrado el estadístico $t$ de la pendiente aterriza en el estadístico $F$ de ANOVA, y los valores críticos también coinciden, $F_{1-\alpha;\,1,n-2} = (t_{1-\alpha/2,\,n-2})^2$. $\blacksquare$

La Figure 10 confirma la mitad distribucional: eleva al cuadrado 200,000 extracciones de una $t_{23}$ y el histograma aterriza sobre la densidad $F_{1,23}$, con el valor crítico del 5% 4.28 igual a 2.069².

El álgebra de arriba ocupa tres líneas, pero la identidad convence mucho más cuando eres tú quien descompone los datos y los dos números siguen negándose a discrepar.

3.8 R cuadrada y su interpretación¶

Intuición¶

La partición ANOVA invita a un resumen natural: ¿qué fracción de la variación total explicó la recta? Esa fracción es $R^2$, el coeficiente de determinación (Definición 3.12). Va de 0 (la recta no explica nada más allá de la media) a 1 (la recta pasa por cada punto). Es el número más reportado y más malinterpretado en la regresión, así que diremos con claridad qué significa y qué no.

Fórmula¶

En palabras: $R^2$ es la proporción de la variación total en $Y$ que la regresión sobre $X$ da cuenta. En la regresión lineal simple también es igual a $r^2$, el cuadrado de la correlación de Pearson entre $X$ y $Y$, un vínculo que el Capítulo 4 desarrolla (4.2 La correlación y la pendiente de regresión).

La Figure 12 muestra los dos ingredientes: la dispersión total de las horas respecto a su media (SSTO, izquierda) se encoge a la dispersión residual mucho menor respecto a la recta (SSE, derecha). $R^2$ es cuánto de la dispersión vertical removió la recta.

3.9 Una banda de confianza para toda la recta¶

Intuición¶

El intervalo de confianza de 3.4 Intervalo de confianza para una respuesta media es honesto sobre un tamaño de lote que nombras de antemano. Pero a menudo quieres trazar una banda alrededor de toda la recta ajustada y afirmar, con 95% de confianza, que la recta verdadera está dentro de ella en todas partes a la vez. Si simplemente encadenas los intervalos puntuales, la confianza simultánea es menor al 95%, porque estás haciendo muchas afirmaciones y algunas fallarán por azar. La banda de Working-Hotelling (Definición 3.13) arregla esto usando un multiplicador un poco mayor, de modo que la afirmación sobre toda la recta se cumpla al nivel enunciado.

Fórmula¶

En palabras: la banda tiene el mismo centro y el mismo error estándar $s\{\hat{Y}_h\}$ que el intervalo puntual, pero el multiplicador $W$ reemplaza a $t_{1-\alpha/2,\,n-2}$. Como $W > t$, la banda es un poco más ancha en todas partes, y ese ancho extra compra cobertura simultánea de la recta verdadera en todos los $X_h$ a la vez. El multiplicador usa una $F$ con 2 grados de libertad en el numerador porque la recta tiene dos parámetros, $\beta_0$ y $\beta_1$, que ambos pueden variar.

W <- sqrt(2 * qf(0.95, 2, n - 2))
round(c(W = W, t = qt(0.975, n - 2)), 4)

     W      t 
2.6162 2.0687 

W = np.sqrt(2 * stats.f.ppf(0.95, 2, n - 2))
print(round(W, 4), round(stats.t.ppf(0.975, n - 2), 4))

2.6162 2.0687

Para Toluca el multiplicador de toda la recta es $W = 2.62$ contra el puntual $t = 2.07$, así que la banda simultánea es alrededor de 26% más ancha que un intervalo de un solo punto. La Figure 13 traza ambos.

3.10 Un ejemplo enfocado en la predicción: publicidad y ventas¶

Para ver la maquinaria en un conjunto de datos nuevo, pasa a los datos de publicidad: los presupuestos de publicidad televisiva de 200 mercados y las ventas de producto resultantes (miles de unidades). Un gerente quiere predecir las ventas de un mercado nuevo que recibirá un presupuesto televisivo de 150 mil dólares, y entender el efecto promedio del presupuesto.

La Figure 14 traza ambas bandas. Nota que la banda de predicción contiene cómodamente casi los 200 mercados, mientras que la banda de confianza, al ser sobre la media, no.

3.11 Resumen del capítulo¶

Ahora puedes hacer inferencia para una regresión lineal simple, no solo ajustarla. Viste por qué la pendiente y el intercepto estandarizados siguen una distribución $t$ con $n-2$ grados de libertad, y la usaste para construir intervalos de confianza y pruebas de hipótesis para ambos coeficientes (la pendiente de Toluca, $t^{\ast} = 10.29$, IC del 95% $(2.85, 4.29)$). Aprendiste a distinguir un intervalo de confianza para una respuesta media de un intervalo de predicción para una sola observación nueva, a derivar ambas fórmulas de varianza, y, sobre todo, a decir cuál necesita una pregunta: para una corrida de 100 unidades el intervalo de la media es $(390, 449)$ pero el intervalo de predicción de una sola corrida es $(314, 525)$. Derivaste la descomposición ANOVA, construiste la tabla, corriste la prueba $F$ ($F^{\ast} = 105.9$), y probaste que es la prueba $t$ de la pendiente al cuadrado. Puedes leer $R^2 = 0.82$ correctamente y ensanchar un intervalo puntual a una banda de Working-Hotelling ($W = 2.62$) cuando necesitas toda la recta.

Resultados clave de un vistazo.

Términos clave. Error estándar estimado, intervalo de confianza para la pendiente, prueba de hipótesis para la pendiente, intervalo de confianza para una respuesta media, intervalo de predicción, análisis de varianza, suma de cuadrados total, suma de cuadrados de la regresión, cuadrado medio, prueba $F$, coeficiente de determinación, banda de confianza de Working-Hotelling, grados de libertad.

Ahora deberías poder:

• Derivar las distribuciones muestrales de $b_0$ y $b_1$ y explicar por qué la inferencia usa $t_{n-2}$.

• Construir e interpretar intervalos de confianza y pruebas de hipótesis para la pendiente y el intercepto a partir de la salida del software.

• Distinguir un intervalo de confianza para una respuesta media de un intervalo de predicción para una observación nueva, derivar ambos y elegir el correcto.

• Derivar la descomposición ANOVA $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$, construir la tabla ANOVA y correr la prueba $F$ para la pendiente.

• Demostrar que $F^{\ast} = (t^{\ast})^2$ en la regresión lineal simple.

• Interpretar $R^2$ correctamente, y nombrar qué mide y qué no mide.

• Construir una banda de confianza de Working-Hotelling para toda la recta de regresión.

Dónde encaja esto. En la columna vertebral del flujo de trabajo de El flujo de trabajo del modelado (PREGUNTAR, EXPLORAR, AJUSTAR, VERIFICAR, USAR), este capítulo es el corazón de la etapa USAR: convierte una recta ajustada en decisiones y pronósticos honestos con incertidumbre declarada. El Capítulo 2 dio las estimaciones puntuales; aquí adjuntamos intervalos y pruebas. El intervalo de predicción (3.5 Intervalo de predicción para una observación nueva) regresa en el Capítulo 12 cuando validamos modelos con datos apartados y en el Capítulo 15 cuando pronosticamos, donde enfrentaremos por qué los datos ordenados en el tiempo hacen más difíciles los intervalos de predicción. La tabla ANOVA (3.6 El enfoque del análisis de varianza) y la prueba $F$ (3.7 La prueba F y su equivalencia con la prueba t) se generalizan en el Capítulo 8 hacia la maquinaria de sumas de cuadrados adicionales (8.3 Sumas de cuadrados adicionales) y la prueba lineal general (8.4 La prueba lineal general), que pueden comparar modelos enteros a la vez. Y el valor p diminuto de la pendiente aquí es la referencia que el Capítulo 5 comprueba contra una prueba de permutación que no supone normalidad en absoluto (5.2 La prueba de permutación para la pendiente), preguntando si nuestra conclusión basada en $t$ sobrevive sin el modelo de errores normales.

3.12 Preguntas frecuentes¶

P1. ¿Cuál es la diferencia en una frase entre un intervalo de confianza y un intervalo de predicción? Un intervalo de confianza es para la respuesta media en $X_h$, promediada sobre todos los casos en ese valor del predictor; un intervalo de predicción es para un solo caso nuevo, y es más ancho porque también debe cubrir la propia dispersión aleatoria de ese caso $\sigma^2$.

P2. ¿Por qué se ensancha el intervalo a medida que $X_h$ se aleja de $\bar{X}$? Porque la varianza de la media ajustada carga un término $(X_h - \bar{X})^2/S_{xx}$. Lejos del centro, un pequeño cambio en la pendiente estimada hace pivotar la recta una gran distancia vertical, así que el valor ajustado ahí es menos certero. Toda recta ajustada pasa por $(\bar{X}, \bar{Y})$, así que la incertidumbre es mínima ahí.

P3. ¿Por qué usar $t$ y no la distribución normal? Porque estimamos $\sigma$ con $s = \sqrt{\mathrm{MSE}}$ en lugar de conocerlo. Esa estimación extra hace que la distribución de referencia sea $t_{n-2}$, con colas más pesadas que la normal. Para $n$ grande las dos casi coinciden, pero para $n = 25$ la diferencia importa.

P4. La prueba $F$ y la prueba $t$ dieron el mismo valor p. ¿Necesito ambas? En la regresión lineal simple, no: $F^{\ast} = (t^{\ast})^2$, así que son una sola prueba. Reporta la prueba $t$ para una sola pendiente porque conserva la dirección (el signo) del efecto. La prueba $F$ se gana su lugar en la regresión múltiple, donde puede probar varios coeficientes de forma simultánea.

P5. Mi $R^2$ es alto. ¿Significa que mi modelo es bueno? No por sí solo. Un $R^2$ alto puede acompañar una relación curva, varianza no constante, o valores atípicos influyentes, todo lo cual un gráfico de residuos revelaría. $R^2$ mide la variación explicada, no la correctitud, ni la exactitud de la predicción, ni la causalidad. Léelo junto con $s$ y los gráficos de diagnóstico (Capítulo 9).

P6. ¿Puedo usar estos intervalos para predecir las horas de un lote de 500 unidades? No. Eso es extrapolación muy fuera del rango de los datos (20 a 120 unidades). Las fórmulas producirán con gusto un intervalo, pero su confianza declarada supone que el modelo de recta se cumple allá afuera, lo cual los datos no pueden respaldar. Confía en los intervalos solo dentro, o apenas fuera, del rango observado de $X$.

P7. ¿Qué significa exactamente “95% de confianza” para el intervalo de la pendiente? Significa que el procedimiento atrapa la pendiente verdadera en el 95% de las muestras repetidas, como muestra la Figure 15 con 100 intervalos simulados. Para tu único intervalo $(2.85, 4.29)$, la pendiente verdadera está en él o no; el 95% describe la confiabilidad a largo plazo del método, no una probabilidad sobre este intervalo en particular.

3.13 Problemas de práctica¶

1. (A) En una frase cada uno, enuncia sobre qué tratan un intervalo de confianza para la respuesta media y un intervalo de predicción para una observación nueva, y di cuál es más ancho y por qué.

2. (A) El intervalo de confianza del 95% de la pendiente es $(2.85, 4.29)$. Explica qué significa “95% de confianza” aquí, y qué no significa.

3. (A) Explica por qué el intervalo de confianza para una respuesta media es más angosto en $X_h = \bar{X}$ y se ensancha a ambos lados.

4. (A) ¿Por qué la inferencia para la pendiente usa una distribución $t$ con $n-2$ grados de libertad en lugar de una normal estándar? ¿De dónde vienen los $n-2$ grados de libertad?

5. (A) Un reporte afirma “$R^2 = 0.82$, así que el modelo es correcto”. Da dos razones distintas por las que esta conclusión no se sigue.

6. (A) El intervalo del 95% del intercepto $(8.21, 116.52)$ es enormemente ancho comparado con el de la pendiente. Explica por qué, refiriéndote a dónde se ubica $X = 0$ respecto a los datos.

7. (A) En la regresión lineal simple la prueba $F$ de ANOVA y la prueba $t$ de la pendiente siempre concuerdan. Enuncia la relación exacta entre los dos estadísticos y una razón para reportar aún la prueba $t$.

8. (A) Un colega usa un intervalo de confianza del 95% para la respuesta media para presupuestar una sola corrida de producción próxima. Explica por qué esto subestima el riesgo, y cuál intervalo debería usar.

9. (B) Deriva el estadístico $t$ $(b_1 - \beta_1)/s\{b_1\} \sim t_{n-2}$ (Teorema 3.2), enunciando con claridad cuáles hechos tomas del Capítulo 2, cuáles tomas por fe hasta el Capítulo 7, y por qué la combinación es una $t$.

10. (B) Deriva $\operatorname{Var}\{\hat{Y}_h\} = \sigma^2(1/n + (X_h - \bar{X})^2/S_{xx})$ (Teorema 3.5), incluyendo el paso de que $\bar{Y}$ y $b_1$ no están correlacionadas.

11. (B) Deriva la varianza del error de predicción $\sigma^2(1 + 1/n + (X_h - \bar{X})^2/S_{xx})$ (Teorema 3.7), explicando de dónde viene el 1 inicial y por qué el error nuevo es independiente de $\hat{Y}_h$.

12. (B) Prueba la identidad ANOVA $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$ (Teorema 3.9) mostrando que el término cruzado se anula, nombrando las identidades de residuos que usas.

13. (B) Prueba que $F^{\ast} = (t^{\ast})^2$ en la regresión lineal simple (Teorema 3.11), partiendo de $\mathrm{SSR} = b_1^2 S_{xx}$.

14. (B) Muestra que $E\{\mathrm{MSR}\} = \sigma^2 + \beta_1^2 S_{xx}$ (Teorema 3.10), y explica por qué esto hace de $\mathrm{MSR}/\mathrm{MSE}$ un estadístico de prueba sensato para $H_0: \beta_1 = 0$.

15. (B) Muestra que un intervalo de confianza del $100(1-\alpha)\%$ para $\beta_1$ excluye a 0 exactamente cuando la prueba $t$ bilateral de $H_0: \beta_1 = 0$ rechaza al nivel $\alpha$.

16. (B) Partiendo de $R^2 = 1 - \mathrm{SSE}/\mathrm{SSTO}$ y $\mathrm{SSE} = S_{yy} - b_1^2 S_{xx}$, muestra que $R^2 = b_1^2 S_{xx}/S_{yy}$, y por lo tanto que $R^2 = r^2$ donde $r = S_{xy}/\sqrt{S_{xx}S_{yy}}$.

17. (B) El multiplicador de Working-Hotelling es $W = \sqrt{2 F_{1-\alpha;2,n-2}}$ y el multiplicador puntual es $t_{1-\alpha/2,n-2}$. Explica, sin una prueba completa, por qué $W$ debe exceder a $t$ y qué compra el ancho extra.

18. (B) Muestra que la semianchura del intervalo de predicción nunca se encoge por debajo de $t\, s$ ni siquiera cuando $n \to \infty$, mientras que la semianchura del intervalo de confianza se encoge a 0. ¿Qué dice esto sobre los límites de la predicción?

19. (C) Ajusta el modelo de Toluca en R o Python. Reproduce el estadístico $t$ de la pendiente, su valor p, y su intervalo de confianza del 95%, y confirma que coinciden con los valores de este capítulo.

20. (C) Calcula un intervalo de confianza del 95% para las horas medias en $X_h = 40$ unidades y en $X_h = 90$ unidades. ¿Cuál es más ancho, y coincide la diferencia con la regla de $(X_h - \bar{X})^2$?

21. (C) Calcula tanto el intervalo de confianza del 95% para la media como el intervalo de predicción del 95% para una corrida nueva en $X_h = 50$ unidades. Reporta el cociente de sus anchos y explícalo.

22. (C) Construye la tabla ANOVA para Toluca en software. De ella, extrae $\mathrm{SSR}$, $\mathrm{SSE}$, $\mathrm{MSR}$, $\mathrm{MSE}$, y $F^{\ast}$, y verifica $F^{\ast} = (t^{\ast})^2$ usando la $t$ de la pendiente del problema 19.

23. (C) Usando los datos de publicidad, ajusta sales ~ TV, prueba la pendiente, y da un intervalo de predicción del 95% para un mercado nuevo con un presupuesto televisivo de 200. Interpreta el intervalo para un gerente.

24. (C) Escribe una simulación con semilla (semilla 4210, 2000 conjuntos de datos del modelo de Toluca ajustado en los tamaños de lote observados) que comprueba la cobertura empírica del intervalo de confianza del 95% para la media en $X_h = 60$. Reporta la cobertura y compárala con 0.95.

25. (C) Reproduce el multiplicador $W$ de Working-Hotelling y el multiplicador puntual $t$ para Toluca. Luego calcula tanto el intervalo puntual del 95% como la semianchura de la banda de Working-Hotelling en $X_h = 120$, y reporta el cociente.

26. (C) Haz la regresión hours ~ lotsize y traza un gráfico de la recta ajustada con la banda de confianza del 95% y la banda de predicción del 95% (como en Figure 1). Describe cómo difieren las dos bandas en ancho y forma.

27. (B) Un estudiante afirma que como $F^{\ast} = 105.9$ es “enorme”, el modelo debe predecir con exactitud las corridas individuales. Explica la confusión entre la prueba $F$ (sobre la pendiente) y la exactitud de la predicción (sobre $s$ y el intervalo de predicción), usando los números de Toluca.

28. (C) Usando los datos de publicidad, calcula $R^2$ de dos maneras: a partir de las sumas de cuadrados de ANOVA y como el cuadrado de la correlación cor(adv$TV, adv$sales). Confirma que coinciden, e interpreta el valor.

3.14 Práctica de examen¶

Los problemas de arriba construyen habilidades una a la vez. Los exámenes te piden combinarlas y, sobre todo, explicar en oraciones completas: leer la salida del software en contexto, juzgar una afirmación, y decir qué cambiaría bajo una condición nueva. Los cinco problemas de abajo están escritos en esa voz de examen. Trabaja cada uno como si fuera para la máxima calificación, en oraciones completas, antes de abrir la respuesta modelo. Un número pelado o un veredicto de una palabra gana poco crédito en el examen real; el razonamiento es la respuesta.

EP 3.1 (interpretar la salida en contexto). Un programador de producción tiene una corrida nueva específica de 80 unidades en el calendario de la próxima semana. El software ajusta el modelo de Toluca y devuelve dos intervalos en un tamaño de lote de 80.

predict(fit, data.frame(lotsize = 80), interval = "confidence")
predict(fit, data.frame(lotsize = 80), interval = "prediction")

      fit      lwr      upr
1 347.982 326.5449 369.4191
      fit      lwr      upr
1 347.982 244.7333 451.2307

pr = fit.get_prediction(pd.DataFrame({"lotsize": [80]}))
print("mean CI ", pr.conf_int(alpha=0.05))
print("pred PI ", pr.conf_int(obs=True, alpha=0.05))

mean CI  [[326.54493825 369.41910215]]
pred PI  [[244.73334785 451.23069256]]

Interpreta cada intervalo en una frase, con unidades y en el contexto de este problema. Luego el programador propone presupuestar la sola corrida de la próxima semana en “entre 327 y 369 horas, con 95% de seguridad”, usando el primer intervalo. Explica en oraciones completas si ese es el intervalo correcto a usar, cuál debería usar el programador en su lugar, y por qué los dos intervalos difieren por tanto.

EP 3.2 (un estudiante afirma X, evalúa). Un compañero de clase escribe: “Bajo el modelo de errores normales la estimación de la pendiente $b_1$ es exactamente normal, así que un intervalo de confianza del 95% para la pendiente debería usar el multiplicador normal estándar 1.96”. El capítulo en cambio usa $t_{0.975,\,23} = 2.069$. Explica en oraciones completas qué está bien y qué está mal en el razonamiento del compañero, qué paso específico fuerza el cambio de 1.96 a 2.069, y qué le pasa a la brecha entre los dos multiplicadores a medida que crece el tamaño de la muestra.

EP 3.3 (interpretar la salida en contexto, evaluar una afirmación). Los datos de publicidad registran los presupuestos televisivos y las ventas de $n = 200$ mercados. Ajustar sales ~ TV da la salida de abajo.

summary(fit_adv)$coefficients
anova(fit_adv)
c(R2 = summary(fit_adv)$r.squared, s = summary(fit_adv)$sigma)

             Estimate Std. Error  t value Pr(>|t|)
(Intercept)  7.032594 0.45784294 15.36028        0
TV           0.047537 0.00269061 17.66763        0
           Df Sum Sq Mean Sq F value    Pr(>F)
TV          1 3314.6  3314.6  312.14 < 2.2e-16
Residuals 198 2102.5    10.6
        R2          s 
    0.6119     3.2587 

print(sm.stats.anova_lm(fit_adv))
print(round(fit_adv.rsquared, 4), round(np.sqrt(fit_adv.scale), 4))

             df       sum_sq      mean_sq           F        PR(>F)
TV          1.0  3314.618167  3314.618167  312.144994  1.467390e-42
Residual  198.0  2102.530583    10.618841         NaN           NaN
0.6119 3.2587

Primero, verifica a partir de los números que el estadístico $F$ es igual al cuadrado del estadístico $t$ de la pendiente de TV, y explica en una frase por qué esa identidad debe cumplirse aquí. Luego un gerente lee $F = 312$ y $R^2 = 0.61$ y concluye: “el modelo es tan significativo que predecirá las ventas de cualquier mercado dentro de una unidad o dos”. Explica en oraciones completas qué establece en realidad el $F$ grande, qué significa $R^2 = 0.61$, y cuál número de la salida gobierna qué tan lejos puede caer de la recta un solo mercado nuevo.

EP 3.4 (qué cambiaría si). Supón que los ingenieros de Toluca reunieran muchos más datos de producción con la misma dispersión de tamaños de lote, de modo que el tamaño de muestra $n$ crece muchísimo mientras $S_{xx}$ crece en proporción y $\mathrm{MSE}$ se mantiene cerca de 2384. Explica en oraciones completas qué le pasa, a medida que $n \to \infty$, a (a) el ancho del intervalo de confianza del 95% para las horas medias en $X_h = 100$, y (b) el ancho del intervalo de predicción del 95% para una sola corrida nueva de 100 unidades. Di por qué los dos se comportan de forma distinta, y qué revela esa diferencia sobre los límites de la predicción.

EP 3.5 (interpretar una figura, explicar por qué). La Figure 5 grafica dos curvas de error estándar estimado contra el tamaño de lote: el error estándar de la respuesta media $s\{\hat{Y}_h\}$ (azul) baja a un mínimo cerca del tamaño de lote medio de 70 y sube hacia ambos extremos, mientras que el error estándar de predicción $s\{\text{pred}\}$ (amarillo) se mantiene casi plano cerca de 49 horas a lo largo de todo el rango y se ubica por encima de la curva azul en todas partes. Explica en oraciones completas (a) por qué la curva azul tiene forma de U con su punto más bajo en $\bar{X}$, y (b) por qué la curva amarilla apenas se dobla aunque está por encima de la azul en cada tamaño de lote.

3.15 Juego del capítulo¶