Estimación de la respuesta, predicción, interpolación y extrapolación

Una vez estimada la recta de regresión, podemos utilizarla para responder preguntas sobre la variable de respuesta \(Y\). Sin embargo, no todas las preguntas se refieren a la misma cantidad.

Debemos distinguir entre:

  1. estimar la respuesta media para un valor determinado de \(X\),
  2. predecir una nueva observación individual para ese mismo valor de \(X\),
  3. realizar la estimación dentro del rango observado, lo que se denomina interpolación,
  4. realizarla fuera del rango observado, lo que se denomina extrapolación.

1. Fundamento conceptual de la predicción

La recta representa una media condicional

El modelo de regresión lineal simple es:

\[ Y_i = \beta_0 + \beta_1X_i + \varepsilon_i \]

La parte sistemática del modelo es:

\[ E(Y\mid X=x) = \beta_0+\beta_1x \]

Después de estimar los parámetros, obtenemos:

\[ \widehat{E(Y\mid X=x)} = \widehat Y = b_0+b_1x \]

La recta de regresión representa la media esperada de \(Y\) para cada valor de \(X\). No representa necesariamente el valor que tomará cada observación individual.

Dos observaciones con el mismo valor de \(X\) pueden presentar valores diferentes de \(Y\), debido al componente aleatorio:

\[ \varepsilon_i \]

En el ejemplo, la recta estimada es:

\[ \widehat{\text{costo}} = 13.3564 + 1.3519 (\text{volumen}) \]

erp_y_ajustada <- fitted(
  erp_modelo
)

plot(
  erp_x,
  erp_y,
  pch = 19,
  cex = 1.15,
  xlab = "Volumen semanal, miles de paquetes",
  ylab = "Costo operativo, miles",
  main = "La recta representa la media condicional"
)

abline(
  erp_modelo,
  lwd = 2.7,
  col = "#1f77b4"
)

segments(
  x0 = erp_x,
  y0 = erp_y,
  x1 = erp_x,
  y1 = erp_y_ajustada,
  col = adjustcolor(
    "#d62728",
    alpha.f = 0.55
  ),
  lwd = 1.5
)

points(
  erp_x,
  erp_y_ajustada,
  pch = 17,
  cex = 1.05,
  col = "#1f77b4"
)

legend(
  "topleft",
  legend = c(
    "Valores observados",
    "Media condicional estimada",
    "Residuales individuales"
  ),
  pch = c(19, 17, NA),
  lty = c(NA, NA, 1),
  lwd = c(NA, NA, 1.5),
  col = c(
    "black",
    "#1f77b4",
    adjustcolor(
      "#d62728",
      alpha.f = 0.65
    )
  ),
  bty = "n"
)

Los triángulos ubicados sobre la recta son las medias estimadas:

\[ \widehat Y_i \]

Los puntos negros son las observaciones individuales:

\[ Y_i \]

La distancia vertical entre ambos es el residual:

\[ e_i=Y_i-\widehat Y_i \]

Distribuciones condicionales alrededor de la recta

La regresión no supone que todas las observaciones se encuentren exactamente sobre la recta. Para cada valor fijo de \(X=x\), existe una distribución de posibles valores de \(Y\) alrededor de la media condicional.

Bajo el modelo normal:

\[ Y\mid X=x \sim N\left( \beta_0+\beta_1x, \sigma^2 \right) \]

Esto significa que:

  • el centro de cada distribución se encuentra sobre la recta,
  • la dispersión vertical representa la variabilidad individual,
  • las observaciones pueden aparecer por encima o por debajo de la media.
erp_x_distribuciones <- as.numeric(
  quantile(
    erp_x,
    probs = c(0.20, 0.50, 0.80),
    names = FALSE
  )
)

erp_y_medias <- predict(
  erp_modelo,
  newdata = data.frame(
    volumen = erp_x_distribuciones
  )
)

erp_limites_y <- range(
  erp_y,
  erp_y_medias - 3.2 * erp_s,
  erp_y_medias + 3.2 * erp_s
)

plot(
  erp_x,
  erp_y,
  pch = 19,
  cex = 1.05,
  xlim = c(
    erp_min_x - 2,
    erp_max_x + 3
  ),
  ylim = erp_limites_y,
  xlab = "Volumen semanal, miles de paquetes",
  ylab = "Costo operativo, miles",
  main = "Distribuciones de Y condicionadas a diferentes valores de X"
)

abline(
  erp_modelo,
  lwd = 2.7,
  col = "#1f77b4"
)

erp_colores_distribuciones <- c(
  "#2e8b57",
  "#7b68ee",
  "#d95f02"
)

for (
  j in seq_along(
    erp_x_distribuciones
  )
) {
  
  erp_y_secuencia <- seq(
    erp_y_medias[j] - 3 * erp_s,
    erp_y_medias[j] + 3 * erp_s,
    length.out = 250
  )
  
  erp_densidad <- dnorm(
    erp_y_secuencia,
    mean = erp_y_medias[j],
    sd = erp_s
  )
  
  erp_densidad_escalada <-
    1.7 *
    erp_densidad /
    max(erp_densidad)
  
  polygon(
    x = c(
      erp_x_distribuciones[j] -
        erp_densidad_escalada,
      rev(
        erp_x_distribuciones[j] +
          erp_densidad_escalada
      )
    ),
    y = c(
      erp_y_secuencia,
      rev(erp_y_secuencia)
    ),
    col = adjustcolor(
      erp_colores_distribuciones[j],
      alpha.f = 0.25
    ),
    border = erp_colores_distribuciones[j],
    lwd = 1.5
  )
  
  points(
    erp_x_distribuciones[j],
    erp_y_medias[j],
    pch = 17,
    cex = 1.35,
    col = erp_colores_distribuciones[j]
  )
  
  segments(
    x0 = erp_x_distribuciones[j] - 1.5,
    y0 = erp_y_medias[j],
    x1 = erp_x_distribuciones[j] + 1.5,
    y1 = erp_y_medias[j],
    lty = 2,
    col = erp_colores_distribuciones[j]
  )
}

legend(
  "topleft",
  legend = c(
    "Observaciones",
    "Media condicional",
    "Distribución de valores individuales"
  ),
  pch = c(19, 17, 15),
  pt.cex = c(1, 1.2, 1.7),
  col = c(
    "black",
    "#1f77b4",
    adjustcolor(
      "#7b68ee",
      alpha.f = 0.45
    )
  ),
  bty = "n"
)

Cada distribución vertical representa posibles valores individuales de \(Y\) para un valor determinado de \(X\).

La recta une los centros de esas distribuciones:

\[ E(Y\mid X=x) \]

Por tanto, debemos distinguir entre el centro de la distribución y una observación obtenida de esa distribución.

2. ¿Qué queremos estimar?

Dos preguntas estadísticas diferentes

Pregunta 1: respuesta media

Supongamos que queremos responder:

¿Cuál es el costo promedio de todas las semanas cuyo volumen es \(x_0\)?

La cantidad objetivo es:

\[ E(Y\mid X=x_0) \]

No estamos tratando de anticipar una semana particular. Estamos estimando el centro de la distribución de costos para todas las semanas con ese volumen.

Para esta pregunta se utiliza un:

\[ \boxed{ \text{intervalo de confianza para la respuesta media} } \]

Pregunta 2: nueva observación individual

Ahora supongamos que queremos responder:

¿Cuál podría ser el costo de una nueva semana específica cuyo volumen sea \(x_0\)?

La cantidad objetivo es:

\[ Y_{\text{nuevo}}\mid X=x_0 \]

En este caso interesa un valor individual que puede quedar por encima o por debajo de la media condicional.

Para esta pregunta se utiliza un:

\[ \boxed{ \text{intervalo de predicción individual} } \]

La predicción puntual es la misma en ambos casos. Lo que cambia es la cantidad que deseamos cubrir y, por tanto, cambia la amplitud del intervalo.

Predicción puntual

Para un valor \(x_0\), la predicción puntual es:

\[ \widehat Y_0 = b_0+b_1x_0 \]

Utilizaremos como ejemplo:

\[ x_0 = 32 \]

Este valor se encuentra dentro del rango observado:

\[ [10, 46] \]

erp_nuevo_inter <- data.frame(
  volumen = erp_x_inter
)

erp_yhat_inter <- as.numeric(
  predict(
    erp_modelo,
    newdata = erp_nuevo_inter
  )
)

cat(
  "Valor de x0:",
  erp_x_inter,
  "\n"
)
#> Valor de x0: 32
cat(
  "Predicción puntual:",
  round(
    erp_yhat_inter,
    5
  ),
  "\n"
)
#> Predicción puntual: 56.618

La predicción puntual:

\[ \widehat Y_0 = 56.618 \]

es simultáneamente:

  • la estimación puntual de la media condicional,
  • el centro del intervalo de confianza,
  • el centro del intervalo de predicción.

3. Construcción de los intervalos

Incertidumbre sobre la media e incertidumbre individual

Factor de posición de \(x_0\)

La incertidumbre depende de la posición de \(x_0\) respecto de la media de \(X\):

\[ h_0 = \frac{1}{n} + \frac{(x_0-\bar X)^2}{S_{XX}} \]

erp_h0_inter <-
  1 / erp_n +
  (
    erp_x_inter -
      erp_x_barra
  )^2 / erp_Sxx

cat(
  "h0 =",
  round(
    erp_h0_inter,
    6
  ),
  "\n"
)
#> h0 = 0.079497

Error estándar de la media estimada

Para estimar la media condicional:

\[ SE_{\text{media}} = s\sqrt{ \frac{1}{n} + \frac{(x_0-\bar X)^2}{S_{XX}} } \]

o, de manera equivalente:

\[ SE_{\text{media}} = s\sqrt{h_0} \]

Error estándar de una nueva observación

Para predecir una observación individual:

\[ SE_{\text{predicción}} = s\sqrt{ 1+ \frac{1}{n} + \frac{(x_0-\bar X)^2}{S_{XX}} } \]

o:

\[ SE_{\text{predicción}} = s\sqrt{1+h_0} \]

erp_SE_media_inter <-
  erp_s *
  sqrt(
    erp_h0_inter
  )

erp_SE_pred_inter <-
  erp_s *
  sqrt(
    1 +
      erp_h0_inter
  )

erp_tabla_errores <- data.frame(
  Cantidad = c(
    "Desviación residual s",
    "Factor h0",
    "Error estándar de la media",
    "Error estándar de predicción"
  ),
  Valor = c(
    erp_s,
    erp_h0_inter,
    erp_SE_media_inter,
    erp_SE_pred_inter
  )
)

knitr::kable(
  erp_tabla_errores,
  digits = 6,
  caption = "Fuentes de incertidumbre para x0"
)
Fuentes de incertidumbre para x0
Cantidad Valor
Desviación residual s 0.636024
Factor h0 0.079497
Error estándar de la media 0.179328
Error estándar de predicción 0.660821

El error estándar de predicción es mayor porque contiene el término adicional:

\[ 1 \]

Este término representa la variabilidad propia de una nueva observación alrededor de su media condicional.

Intervalo de confianza para la respuesta media

El intervalo de confianza de nivel \(1-\alpha\) es:

\[ \widehat Y_0 \pm t_{\alpha/2,n-2} s \sqrt{ \frac{1}{n} + \frac{(x_0-\bar X)^2}{S_{XX}} } \]

erp_IC_media_inter <- c(
  inferior =
    erp_yhat_inter -
    erp_t_critico *
    erp_SE_media_inter,
  
  superior =
    erp_yhat_inter +
    erp_t_critico *
    erp_SE_media_inter
)

erp_IC_media_inter
#> inferior superior 
#>   56.227   57.009

Este intervalo responde:

¿En qué rango se encuentra el costo medio de todas las semanas cuyo volumen es \(x_0\)?

Intervalo de predicción individual

El intervalo de predicción de nivel \(1-\alpha\) es:

\[ \widehat Y_0 \pm t_{\alpha/2,n-2} s \sqrt{ 1+ \frac{1}{n} + \frac{(x_0-\bar X)^2}{S_{XX}} } \]

erp_IP_individual_inter <- c(
  inferior =
    erp_yhat_inter -
    erp_t_critico *
    erp_SE_pred_inter,
  
  superior =
    erp_yhat_inter +
    erp_t_critico *
    erp_SE_pred_inter
)

erp_IP_individual_inter
#> inferior superior 
#>   55.178   58.058

Este intervalo responde:

¿En qué rango podría encontrarse el costo de una nueva semana específica cuyo volumen es \(x_0\)?

Comparación numérica

erp_tabla_intervalos_inter <- data.frame(
  Intervalo = c(
    "Confianza para la media",
    "Predicción individual"
  ),
  Objetivo = c(
    "E(Y | X = x0)",
    "Y nuevo | X = x0"
  ),
  Estimacion = c(
    erp_yhat_inter,
    erp_yhat_inter
  ),
  Limite_inferior = c(
    erp_IC_media_inter["inferior"],
    erp_IP_individual_inter["inferior"]
  ),
  Limite_superior = c(
    erp_IC_media_inter["superior"],
    erp_IP_individual_inter["superior"]
  )
)

erp_tabla_intervalos_inter$Amplitud <-
  erp_tabla_intervalos_inter$Limite_superior -
  erp_tabla_intervalos_inter$Limite_inferior

knitr::kable(
  erp_tabla_intervalos_inter,
  digits = 5,
  caption = paste(
    "Intervalos para x0 =",
    round(
      erp_x_inter,
      2
    )
  )
)
Intervalos para x0 = 32
Intervalo Objetivo Estimacion Limite_inferior Limite_superior Amplitud
Confianza para la media E(Y | X = x0) 56.618 56.227 57.009 0.78144
Predicción individual Y nuevo | X = x0 56.618 55.178 58.058 2.87961

Verificación mediante predict()

erp_resultado_confianza <- predict(
  erp_modelo,
  newdata = erp_nuevo_inter,
  interval = "confidence",
  level = erp_nivel
)

erp_resultado_prediccion <- predict(
  erp_modelo,
  newdata = erp_nuevo_inter,
  interval = "prediction",
  level = erp_nivel
)

erp_resultado_confianza
#>      fit    lwr    upr
#> 1 56.618 56.227 57.009
erp_resultado_prediccion
#>      fit    lwr    upr
#> 1 56.618 55.178 58.058

4. Interpretación gráfica de la incertidumbre

Distribución de la media estimada y distribución predictiva

Las dos distribuciones están centradas en la misma predicción puntual:

\[ \widehat Y_0 \]

pero poseen dispersiones diferentes.

La distribución muestral aproximada de la media estimada utiliza:

\[ SE_{\text{media}} \]

La distribución predictiva aproximada de una nueva observación utiliza:

\[ SE_{\text{predicción}} \]

Comparación directa

erp_secuencia_y <- seq(
  erp_yhat_inter -
    4.2 *
    erp_SE_pred_inter,
  
  erp_yhat_inter +
    4.2 *
    erp_SE_pred_inter,
  
  length.out = 600
)

erp_densidad_media <- dnorm(
  erp_secuencia_y,
  mean = erp_yhat_inter,
  sd = erp_SE_media_inter
)

erp_densidad_pred <- dnorm(
  erp_secuencia_y,
  mean = erp_yhat_inter,
  sd = erp_SE_pred_inter
)

plot(
  erp_secuencia_y,
  erp_densidad_pred,
  type = "l",
  lwd = 2.8,
  col = "#d95f02",
  xlab = "Posibles valores del costo",
  ylab = "Densidad",
  main = paste(
    "Confianza y predicción para x0 =",
    round(
      erp_x_inter,
      2
    )
  )
)

lines(
  erp_secuencia_y,
  erp_densidad_media,
  lwd = 2.8,
  col = "#1f77b4"
)

abline(
  v = erp_yhat_inter,
  lty = 2,
  lwd = 2,
  col = "#444444"
)

abline(
  v = erp_IC_media_inter,
  lty = 3,
  lwd = 1.7,
  col = "#1f77b4"
)

abline(
  v = erp_IP_individual_inter,
  lty = 3,
  lwd = 1.7,
  col = "#d95f02"
)

legend(
  "topright",
  legend = c(
    "Distribución de la media estimada",
    "Distribución de una observación nueva",
    "Predicción puntual"
  ),
  lty = c(1, 1, 2),
  lwd = c(2.8, 2.8, 2),
  col = c(
    "#1f77b4",
    "#d95f02",
    "#444444"
  ),
  bty = "n"
)

La curva azul es más estrecha porque únicamente representa la incertidumbre al estimar el centro de la distribución.

La curva naranja es más ancha porque representa dónde podría aparecer una observación individual.

Áreas centrales de 95%

par(
  mfrow = c(1, 2)
)

# Intervalo de confianza

plot(
  erp_secuencia_y,
  erp_densidad_media,
  type = "l",
  lwd = 2.5,
  col = "#1f77b4",
  xlab = "Costo medio posible",
  ylab = "Densidad",
  main = "IC de 95% para la media"
)

erp_indice_ic <-
  erp_secuencia_y >=
  erp_IC_media_inter["inferior"] &
  erp_secuencia_y <=
  erp_IC_media_inter["superior"]

polygon(
  x = c(
    erp_IC_media_inter["inferior"],
    erp_secuencia_y[erp_indice_ic],
    erp_IC_media_inter["superior"]
  ),
  y = c(
    0,
    erp_densidad_media[erp_indice_ic],
    0
  ),
  col = adjustcolor(
    "#1f77b4",
    alpha.f = 0.35
  ),
  border = NA
)

lines(
  erp_secuencia_y,
  erp_densidad_media,
  lwd = 2.5,
  col = "#1f77b4"
)

abline(
  v = erp_yhat_inter,
  lty = 2
)

# Intervalo de predicción

plot(
  erp_secuencia_y,
  erp_densidad_pred,
  type = "l",
  lwd = 2.5,
  col = "#d95f02",
  xlab = "Costo individual posible",
  ylab = "Densidad",
  main = "IP de 95% para un individuo"
)

erp_indice_ip <-
  erp_secuencia_y >=
  erp_IP_individual_inter["inferior"] &
  erp_secuencia_y <=
  erp_IP_individual_inter["superior"]

polygon(
  x = c(
    erp_IP_individual_inter["inferior"],
    erp_secuencia_y[erp_indice_ip],
    erp_IP_individual_inter["superior"]
  ),
  y = c(
    0,
    erp_densidad_pred[erp_indice_ip],
    0
  ),
  col = adjustcolor(
    "#d95f02",
    alpha.f = 0.35
  ),
  border = NA
)

lines(
  erp_secuencia_y,
  erp_densidad_pred,
  lwd = 2.5,
  col = "#d95f02"
)

abline(
  v = erp_yhat_inter,
  lty = 2
)

par(
  mfrow = c(1, 1)
)

En una interpretación frecuentista, el 95% se refiere al desempeño del procedimiento en repeticiones del muestreo. No significa que, después de calcular los límites, exista 95% de probabilidad de que un parámetro fijo esté dentro del intervalo.

Bandas de confianza y predicción sobre la recta

Cuando los intervalos se calculan para muchos valores de \(X\), obtenemos bandas alrededor de la recta.

erp_grid_dentro <- data.frame(
  volumen = seq(
    erp_min_x,
    erp_max_x,
    length.out = 300
  )
)

erp_banda_confianza <- predict(
  erp_modelo,
  newdata = erp_grid_dentro,
  interval = "confidence",
  level = erp_nivel
)

erp_banda_prediccion <- predict(
  erp_modelo,
  newdata = erp_grid_dentro,
  interval = "prediction",
  level = erp_nivel
)

erp_limites_bandas <- range(
  erp_y,
  erp_banda_prediccion[
    ,
    "lwr"
  ],
  erp_banda_prediccion[
    ,
    "upr"
  ]
)

plot(
  erp_x,
  erp_y,
  pch = 19,
  cex = 1.05,
  xlim = range(
    erp_grid_dentro$volumen
  ),
  ylim = erp_limites_bandas,
  xlab = "Volumen semanal, miles de paquetes",
  ylab = "Costo operativo, miles",
  main = "Bandas de confianza y predicción"
)

polygon(
  x = c(
    erp_grid_dentro$volumen,
    rev(
      erp_grid_dentro$volumen
    )
  ),
  y = c(
    erp_banda_prediccion[
      ,
      "lwr"
    ],
    rev(
      erp_banda_prediccion[
        ,
        "upr"
      ]
    )
  ),
  col = adjustcolor(
    "#f4a261",
    alpha.f = 0.28
  ),
  border = NA
)

polygon(
  x = c(
    erp_grid_dentro$volumen,
    rev(
      erp_grid_dentro$volumen
    )
  ),
  y = c(
    erp_banda_confianza[
      ,
      "lwr"
    ],
    rev(
      erp_banda_confianza[
        ,
        "upr"
      ]
    )
  ),
  col = adjustcolor(
    "#4c78a8",
    alpha.f = 0.38
  ),
  border = NA
)

lines(
  erp_grid_dentro$volumen,
  erp_banda_confianza[
    ,
    "fit"
  ],
  lwd = 2.7,
  col = "#1f4e79"
)

points(
  erp_x,
  erp_y,
  pch = 19
)

points(
  erp_x_inter,
  erp_yhat_inter,
  pch = 17,
  cex = 1.5,
  col = "#b22222"
)

legend(
  "topleft",
  legend = c(
    "Recta estimada",
    "IC para la media",
    "IP individual",
    "Predicción puntual"
  ),
  lty = c(
    1,
    NA,
    NA,
    NA
  ),
  lwd = c(
    2.7,
    NA,
    NA,
    NA
  ),
  pch = c(
    NA,
    15,
    15,
    17
  ),
  pt.cex = c(
    NA,
    2,
    2,
    1.3
  ),
  col = c(
    "#1f4e79",
    adjustcolor(
      "#4c78a8",
      alpha.f = 0.65
    ),
    adjustcolor(
      "#f4a261",
      alpha.f = 0.65
    ),
    "#b22222"
  ),
  bty = "n"
)

La banda azul representa intervalos para la media condicional.

La banda naranja representa intervalos para observaciones individuales.

La banda de predicción es más amplia en todos los valores de \(X\).

¿Por qué los intervalos se ensanchan?

La expresión:

\[ \frac{(x_0-\bar X)^2}{S_{XX}} \]

mide la distancia de \(x_0\) respecto del centro de los datos.

Cuando:

\[ x_0=\bar X \]

esa cantidad es igual a cero y los intervalos alcanzan su menor amplitud.

Al alejarse de \(\bar X\), la incertidumbre aumenta.

erp_ancho_ic <-
  erp_banda_confianza[
    ,
    "upr"
  ] -
  erp_banda_confianza[
    ,
    "lwr"
  ]

erp_ancho_ip <-
  erp_banda_prediccion[
    ,
    "upr"
  ] -
  erp_banda_prediccion[
    ,
    "lwr"
  ]

matplot(
  erp_grid_dentro$volumen,
  cbind(
    erp_ancho_ic,
    erp_ancho_ip
  ),
  type = "l",
  lty = c(1, 2),
  lwd = 2.7,
  col = c(
    "#1f77b4",
    "#d95f02"
  ),
  xlab = "Valor de X",
  ylab = "Amplitud total del intervalo",
  main = "Amplitud de los intervalos según la posición de X"
)

abline(
  v = erp_x_barra,
  lty = 3,
  lwd = 2,
  col = "#444444"
)

legend(
  "top",
  legend = c(
    "IC para la media",
    "IP individual",
    "Media de X"
  ),
  lty = c(1, 2, 3),
  lwd = c(2.7, 2.7, 2),
  col = c(
    "#1f77b4",
    "#d95f02",
    "#444444"
  ),
  bty = "n"
)

5. Interpolación y extrapolación

Interpolación

Se denomina interpolación al uso del modelo para un valor de \(X\) que se encuentra dentro del rango observado.

Formalmente:

\[ \min(X) \leq x_0 \leq \max(X) \]

En este ejemplo, el rango observado es:

\[ 10 \leq X \leq 46 \]

Como:

\[ x_0 = 32 \]

se encuentra dentro de ese rango, se trata de una interpolación.

La interpolación utiliza información ubicada entre valores que efectivamente fueron observados. Por esa razón, suele ser más defendible que la extrapolación.

Extrapolación

Se denomina extrapolación al uso del modelo para un valor de \(X\) que se encuentra fuera del rango observado.

Formalmente:

\[ x_0<\min(X) \]

o:

\[ x_0>\max(X) \]

En este ejemplo utilizaremos:

\[ x_{\text{extra}} = 53.2 \]

Como este valor es mayor que:

\[ \max(X) = 46 \]

la estimación constituye una extrapolación.

R puede calcular predicciones e intervalos fuera del rango observado. Sin embargo, obtener un resultado numérico no garantiza que la relación lineal continúe siendo válida en esa región.

Gráfico de interpolación y extrapolación

erp_extension_izquierda <-
  0.10 *
  erp_rango_x

erp_extension_derecha <-
  max(
    erp_x_extra -
      erp_max_x,
    0.20 *
      erp_rango_x
  )

erp_grid_ampliado <- data.frame(
  volumen = seq(
    erp_min_x -
      erp_extension_izquierda,
    erp_max_x +
      erp_extension_derecha,
    length.out = 400
  )
)

erp_conf_ampliado <- predict(
  erp_modelo,
  newdata = erp_grid_ampliado,
  interval = "confidence",
  level = erp_nivel
)

erp_pred_ampliado <- predict(
  erp_modelo,
  newdata = erp_grid_ampliado,
  interval = "prediction",
  level = erp_nivel
)

erp_limites_y_ampliados <- range(
  erp_y,
  erp_pred_ampliado[
    ,
    "lwr"
  ],
  erp_pred_ampliado[
    ,
    "upr"
  ]
)

plot(
  erp_x,
  erp_y,
  pch = 19,
  cex = 1.05,
  xlim = range(
    erp_grid_ampliado$volumen
  ),
  ylim = erp_limites_y_ampliados,
  xlab = "Volumen semanal, miles de paquetes",
  ylab = "Costo operativo, miles",
  main = "Interpolación y extrapolación"
)

rect(
  xleft = erp_min_x,
  ybottom = par("usr")[3],
  xright = erp_max_x,
  ytop = par("usr")[4],
  col = adjustcolor(
    "#90c987",
    alpha.f = 0.12
  ),
  border = NA
)

polygon(
  x = c(
    erp_grid_ampliado$volumen,
    rev(
      erp_grid_ampliado$volumen
    )
  ),
  y = c(
    erp_pred_ampliado[
      ,
      "lwr"
    ],
    rev(
      erp_pred_ampliado[
        ,
        "upr"
      ]
    )
  ),
  col = adjustcolor(
    "#f4a261",
    alpha.f = 0.23
  ),
  border = NA
)

polygon(
  x = c(
    erp_grid_ampliado$volumen,
    rev(
      erp_grid_ampliado$volumen
    )
  ),
  y = c(
    erp_conf_ampliado[
      ,
      "lwr"
    ],
    rev(
      erp_conf_ampliado[
        ,
        "upr"
      ]
    )
  ),
  col = adjustcolor(
    "#4c78a8",
    alpha.f = 0.32
  ),
  border = NA
)

lines(
  erp_grid_ampliado$volumen,
  erp_conf_ampliado[
    ,
    "fit"
  ],
  lwd = 2.7,
  col = "#1f4e79"
)

points(
  erp_x,
  erp_y,
  pch = 19
)

abline(
  v = erp_min_x,
  lty = 2,
  lwd = 2,
  col = "#555555"
)

abline(
  v = erp_max_x,
  lty = 2,
  lwd = 2,
  col = "#555555"
)

abline(
  v = erp_x_inter,
  lty = 3,
  lwd = 2,
  col = "#2e8b57"
)

abline(
  v = erp_x_extra,
  lty = 3,
  lwd = 2,
  col = "#b22222"
)

erp_y_inter_graf <- predict(
  erp_modelo,
  newdata = data.frame(
    volumen = erp_x_inter
  )
)

erp_y_extra_graf <- predict(
  erp_modelo,
  newdata = data.frame(
    volumen = erp_x_extra
  )
)

points(
  erp_x_inter,
  erp_y_inter_graf,
  pch = 17,
  cex = 1.5,
  col = "#2e8b57"
)

points(
  erp_x_extra,
  erp_y_extra_graf,
  pch = 17,
  cex = 1.5,
  col = "#b22222"
)

text(
  x = (
    erp_min_x +
      erp_max_x
  ) / 2,
  y = par("usr")[4] -
    0.06 *
    diff(
      par("usr")[3:4]
    ),
  labels = "Zona de interpolación",
  col = "#2e6b31",
  font = 2
)

text(
  x = erp_x_extra,
  y = erp_y_extra_graf,
  labels = "  extrapolación",
  pos = 4,
  col = "#b22222"
)

legend(
  "topleft",
  legend = c(
    "Rango observado",
    "Recta estimada",
    "IC para la media",
    "IP individual",
    "Interpolación",
    "Extrapolación"
  ),
  pch = c(
    15,
    NA,
    15,
    15,
    17,
    17
  ),
  lty = c(
    NA,
    1,
    NA,
    NA,
    NA,
    NA
  ),
  lwd = c(
    NA,
    2.7,
    NA,
    NA,
    NA,
    NA
  ),
  pt.cex = c(
    2,
    NA,
    2,
    2,
    1.2,
    1.2
  ),
  col = c(
    adjustcolor(
      "#90c987",
      alpha.f = 0.65
    ),
    "#1f4e79",
    adjustcolor(
      "#4c78a8",
      alpha.f = 0.65
    ),
    adjustcolor(
      "#f4a261",
      alpha.f = 0.65
    ),
    "#2e8b57",
    "#b22222"
  ),
  bty = "n"
)

¿Qué intervalo debe utilizarse?

El tipo de intervalo depende de la cantidad objetivo, no solamente de si estamos interpolando o extrapolando.

Si queremos una media

Para responder:

¿Cuál es la respuesta promedio esperada cuando \(X=x_0\)?

utilizamos:

\[ \boxed{ \text{intervalo de confianza para la media} } \]

Esto puede calcularse tanto en interpolación como en extrapolación, aunque en extrapolación debe interpretarse con cautela.

Si queremos una observación individual

Para responder:

¿Cuál podría ser el valor de una nueva observación cuando \(X=x_0\)?

utilizamos:

\[ \boxed{ \text{intervalo de predicción individual} } \]

También puede calcularse dentro o fuera del rango observado, pero la extrapolación implica un riesgo adicional sobre la validez del modelo.

Tabla de decisión

erp_tabla_decision <- data.frame(
  Posicion_de_x0 = c(
    "Dentro del rango observado",
    "Dentro del rango observado",
    "Fuera del rango observado",
    "Fuera del rango observado"
  ),
  Situacion = c(
    "Interpolación",
    "Interpolación",
    "Extrapolación",
    "Extrapolación"
  ),
  Pregunta = c(
    "Media esperada de Y",
    "Nueva observación individual",
    "Media esperada de Y",
    "Nueva observación individual"
  ),
  Intervalo = c(
    "Intervalo de confianza",
    "Intervalo de predicción",
    "Intervalo de confianza",
    "Intervalo de predicción"
  ),
  Recomendacion = c(
    "Uso habitual",
    "Uso habitual",
    "Interpretar con cautela",
    "Interpretar con mucha cautela"
  )
)

knitr::kable(
  erp_tabla_decision,
  caption = "Selección del intervalo según la pregunta y la posición de x0"
)
Selección del intervalo según la pregunta y la posición de x0
Posicion_de_x0 Situacion Pregunta Intervalo Recomendacion
Dentro del rango observado Interpolación Media esperada de Y Intervalo de confianza Uso habitual
Dentro del rango observado Interpolación Nueva observación individual Intervalo de predicción Uso habitual
Fuera del rango observado Extrapolación Media esperada de Y Intervalo de confianza Interpretar con cautela
Fuera del rango observado Extrapolación Nueva observación individual Intervalo de predicción Interpretar con mucha cautela

Interpolación y extrapolación describen la posición de \(x_0\). Intervalo de confianza e intervalo de predicción describen la cantidad que queremos estimar. Son clasificaciones distintas.

6. Ejemplo aplicado y reglas de decisión

Ejemplo consolidado

Compararemos dos valores:

\[ x_{\text{inter}} = 32 \]

que se encuentra dentro del rango observado, y:

\[ x_{\text{extra}} = 53.2 \]

que se encuentra fuera del rango observado.

erp_nuevos_ejemplo <- data.frame(
  volumen = c(
    erp_x_inter,
    erp_x_extra
  )
)

erp_conf_ejemplo <- predict(
  erp_modelo,
  newdata = erp_nuevos_ejemplo,
  interval = "confidence",
  level = erp_nivel
)

erp_pred_ejemplo <- predict(
  erp_modelo,
  newdata = erp_nuevos_ejemplo,
  interval = "prediction",
  level = erp_nivel
)

erp_clasificar_posicion <- function(
  valor_x,
  minimo_x,
  maximo_x
) {
  ifelse(
    valor_x >= minimo_x &
      valor_x <= maximo_x,
    "Interpolación",
    "Extrapolación"
  )
}

erp_tabla_ejemplo <- data.frame(
  volumen = erp_nuevos_ejemplo$volumen,
  situacion = erp_clasificar_posicion(
    erp_nuevos_ejemplo$volumen,
    erp_min_x,
    erp_max_x
  ),
  prediccion_puntual = erp_conf_ejemplo[
    ,
    "fit"
  ],
  IC_media_inferior = erp_conf_ejemplo[
    ,
    "lwr"
  ],
  IC_media_superior = erp_conf_ejemplo[
    ,
    "upr"
  ],
  IP_individual_inferior = erp_pred_ejemplo[
    ,
    "lwr"
  ],
  IP_individual_superior = erp_pred_ejemplo[
    ,
    "upr"
  ]
)

knitr::kable(
  erp_tabla_ejemplo,
  digits = 4,
  caption = "Intervalos en interpolación y extrapolación"
)
Intervalos en interpolación y extrapolación
volumen situacion prediccion_puntual IC_media_inferior IC_media_superior IP_individual_inferior IP_individual_superior
32.0 Interpolación 56.618 56.227 57.009 55.178 58.058
53.2 Extrapolación 85.279 84.378 86.180 83.626 86.932

Interpretación de la interpolación

Para:

\[ X = 32 \]

la predicción puntual es:

\[ \widehat Y = 56.618 \]

El intervalo de confianza para la media es:

\[ [ 56.2273, 57.0087 ] \]

Este intervalo se utiliza para estimar el costo medio de todas las semanas con ese volumen.

El intervalo de predicción individual es:

\[ [ 55.1782, 58.0578 ] \]

Este intervalo se utiliza para anticipar el costo de una nueva semana específica con ese volumen.

Interpretación de la extrapolación

Para:

\[ X = 53.2 \]

la predicción puntual es:

\[ \widehat Y = 85.2788 \]

El intervalo de confianza para la media es:

\[ [ 84.3776, 86.1799 ] \]

El intervalo de predicción individual es:

\[ [ 83.6258, 86.9318 ] \]

Los cálculos son matemáticamente válidos bajo el modelo, pero su interpretación requiere cautela porque el valor de \(X\) está fuera del rango observado.

¿Por qué la extrapolación es riesgosa?

Cuando extrapolamos, asumimos que fuera del rango observado:

  1. la relación continúa siendo lineal,
  2. la pendiente continúa siendo la misma,
  3. la variabilidad residual no cambia,
  4. no aparecen límites físicos, financieros u operativos,
  5. no existen variables omitidas que cambien la relación.

El intervalo aumenta al alejarse de \(\bar X\), pero ese ensanchamiento únicamente refleja la incertidumbre calculada bajo el mismo modelo lineal.

No incorpora automáticamente la posibilidad de que la forma funcional sea incorrecta fuera del rango observado.

Un intervalo más amplio no elimina el riesgo de extrapolación. El modelo sigue suponiendo que la recta continúa siendo válida fuera de la región donde se obtuvieron los datos.

¿Puede el intervalo salir del rango observado de Y?

Sí. Esto puede ocurrir incluso durante una interpolación.

Supongamos que los valores observados de \(Y\) se encuentran entre:

\[ \min(Y) = 27 \]

y:

\[ \max(Y) = 76 \]

Un intervalo de predicción puede producir límites menores que el mínimo observado o mayores que el máximo observado.

Esto no representa necesariamente un error, porque el intervalo se refiere a posibles observaciones futuras, no solamente a los valores que ya aparecieron en la muestra.

erp_comparacion_rangos <- data.frame(
  Cantidad = c(
    "Mínimo observado de Y",
    "Máximo observado de Y",
    "Límite inferior del IP en interpolación",
    "Límite superior del IP en interpolación",
    "Límite inferior del IP en extrapolación",
    "Límite superior del IP en extrapolación"
  ),
  Valor = c(
    min(erp_y),
    max(erp_y),
    erp_pred_ejemplo[1, "lwr"],
    erp_pred_ejemplo[1, "upr"],
    erp_pred_ejemplo[2, "lwr"],
    erp_pred_ejemplo[2, "upr"]
  )
)

knitr::kable(
  erp_comparacion_rangos,
  digits = 4,
  caption = "Comparación entre el rango observado de Y y los intervalos"
)
Comparación entre el rango observado de Y y los intervalos
Cantidad Valor
Mínimo observado de Y 27.000
Máximo observado de Y 76.000
Límite inferior del IP en interpolación 55.178
Límite superior del IP en interpolación 58.058
Límite inferior del IP en extrapolación 83.626
Límite superior del IP en extrapolación 86.932

Preguntas para distinguir los conceptos

Situación A

¿Cuál es el costo promedio de todas las semanas que procesan 32 mil paquetes?

La pregunta se refiere a una media:

\[ E(Y\mid X=32) \]

Como 32 está dentro del rango observado:

  • situación: interpolación,
  • intervalo: intervalo de confianza para la media.

Situación B

¿Cuál podría ser el costo de la próxima semana si se procesan 32 mil paquetes?

La pregunta se refiere a una observación individual:

\[ Y_{\text{nuevo}}\mid X=32 \]

Como 32 está dentro del rango observado:

  • situación: interpolación,
  • intervalo: intervalo de predicción individual.

Situación C

¿Cuál sería el costo promedio si el volumen aumentara a un nivel superior al máximo registrado?

La pregunta se refiere a una media, pero el valor de \(X\) está fuera del rango:

  • situación: extrapolación,
  • intervalo: intervalo de confianza para la media,
  • advertencia: la interpretación depende de que la recta siga siendo válida.

Situación D

¿Cuál podría ser el costo de una nueva semana con un volumen superior al máximo registrado?

La pregunta se refiere a una observación individual y el valor de \(X\) está fuera del rango:

  • situación: extrapolación,
  • intervalo: intervalo de predicción individual,
  • advertencia: es la situación que requiere mayor cautela.

Regla práctica final

La selección puede realizarse mediante dos preguntas sucesivas.

Primera pregunta: ¿dónde está \(x_0\)?

Si:

\[ \min(X) \leq x_0 \leq \max(X) \]

se trata de:

\[ \boxed{\text{interpolación}} \]

Si:

\[ x_0<\min(X) \]

o:

\[ x_0>\max(X) \]

se trata de:

\[ \boxed{\text{extrapolación}} \]

Segunda pregunta: ¿qué queremos estimar?

Si queremos:

\[ E(Y\mid X=x_0) \]

utilizamos:

\[ \boxed{\text{intervalo de confianza}} \]

Si queremos:

\[ Y_{\text{nuevo}}\mid X=x_0 \]

utilizamos:

\[ \boxed{\text{intervalo de predicción}} \]

Síntesis

\[ \text{posición de }x_0 \longrightarrow \begin{cases} \text{interpolación}\\ \text{extrapolación} \end{cases} \]

\[ \text{cantidad objetivo} \longrightarrow \begin{cases} \text{media} \rightarrow \text{intervalo de confianza}\\ \text{individuo} \rightarrow \text{intervalo de predicción} \end{cases} \]

La interpolación o extrapolación indica dónde se encuentra \(x_0\). El intervalo de confianza o de predicción indica qué cantidad se desea cubrir. Una clasificación no sustituye a la otra.