Intervalos de predicción

El intervalo de predicción para \(\sum_{j=1}^{m}a_{j}y_{hj}\) es

\[ \left[ \hat{\beta_{0}}\sum_{j=1}^{n}a_{j}+\hat{\beta_{1}}\sum_{j=1}^{n}a_{j}x_{hj}\pm t_{n-2,1-\alpha/2 }\sqrt{\hat{\sigma}^{2}\left[\sum_{j=1}^{m}a_{j}^{2}+\frac{\left(\sum_{j=1}^{m} a_{j}\right)^{2}}{n} +\frac{\left(\sum_{j=1}^{m}a_{j}x_{hj}-\bar{X}\sum_{j=1}^{m}a_{j}\right)}{Ssx}\right]} \right] \]

Continuaremos con los datos de la compañía Toluca. Además, calcularemos algunas estadísticas para el intervalo de predicción, donde recordemos que no está automatizado en R para el caso de las combinaciones lineales de los \(y_{h_{j}}\)

library(ALSM)
library(latex2exp)
Datos=TolucaCompany

# Cálculos: 
xbar=mean(Datos$x)
SSx=sum((Datos$x-xbar)^2)
ybar=mean(Datos$y)
SSy=sum((Datos$y-ybar)^2)
SSxy=sum((Datos$y-ybar)*(Datos$x-xbar))
beta1=SSxy/SSx
beta0=ybar-beta1*xbar
n=length(Datos$x)
Datos$yhat=beta0+beta1*Datos$x
Datos$error=Datos$y-Datos$yhat
MSE=sum((Datos$error)^2)/(n-2)

y ajustamos el modelo, el cual ya sabemos que todos los supuestos se cumplen

fit=lm(y~x, data=Datos)

Ahora bien:

\(a)\) Suponga que el siguiente lote que se producirá será de tamaño 100, el directivo está interesado en saber cuántas horas se trabajarán para la producción de ese lote. Dé un intervalo al 90% de predicción.

Solución: Cuando sólo se realiza la predicción de una sóla observación se puede usar directamente la función predict en R con interval = "prediction". Para ello

# La observación a predecir la guardamos en un dataframe, donde x corresponde al nombre de la 
# columna de las observaciones
newdata <- data.frame(x = c(100))

# Realizamos la predicción al 90% de predicción 
Predx100 <- predict(fit, newdata, interval = "prediction", level = 0.90)
head(Predx100)
##        fit      lwr      upr
## 1 419.3861 332.2072 506.5649

siendo fit la estimación puntual de la predicción y lwr y upr los límites inferior y superior, respectivamente, del intervalo de predicción. De tal manera podemos responder a la pregunta diciendo que que las horas que se requieren para trabajar en este nuevo lote están entre las 332.2 y 506.56 horas con una confianza del 90%.

Alternativa

Podemos realizar los “cálculos a mano” siguiendo las fórmula de las notas

var.PredYh_x100=MSE*(1+1/n+(100-xbar)^2/SSx)
alpha=.1
tcuantil=qt(1-alpha/2, n-2,lower.tail = TRUE)
print(c(beta0+beta1*100-tcuantil*sqrt(var.PredYh_x100),
        beta0+beta1*100+tcuantil*sqrt(var.PredYh_x100)))
## [1] 332.2072 506.5649

que nos arroja el mismo intervalo de predicción.

Observación: Notar que los intervalos de predicción son en general más amplios que los asociados a la \(\mathbb{E}[y|x]\):

predict(fit, newdata, interval = "confidence", level = 0.90)
##        fit      lwr     upr
## 1 419.3861 394.9251 443.847

\(b)\) Suponga que la compañía ha firmado un contrato para producir en tres de sus unidades, un lote de tamaño 100 en cada uno. El directivo está interesado en saber cuántas horas en promedio por unidad se usarán para cumplir con lo establecido en el contrato.

Solución: En este caso sólo se puede a mano, buscamos la estimación de \(\sum_{j=1}^3 \frac{y_{hj}}{3}\), con \(x_{h1}= x_{h2}=x_{h3}=100\). De donde \(a_j=1/3\) para \(j=1,2,3\). Realizamos a “mano” los cálculos

sumA_j=(1/3+1/3+1/3)
sumA_j_2=(1/3^2+1/3^2+1/3^2)
SumA_jXhj=((1/3)*100+(1/3)*100+(1/3)*100)
var.PredMeanYh_x100=MSE*(sumA_j_2+sumA_j^2/n+(SumA_jXhj-sumA_j*xbar)^2/SSx)
alpha=.1
tcuantil=qt(1-alpha/2, n-2,lower.tail = TRUE)
print(c(sumA_j*beta0+beta1*SumA_jXhj-tcuantil*sqrt(var.PredMeanYh_x100),
        sumA_j*beta0+beta1*SumA_jXhj+tcuantil*sqrt(var.PredMeanYh_x100)))
## [1] 365.2356 473.5366

En promedio se van a requerir para estos tres lotes, por cada uno, entre 365.2 y 473.5 horas Notar que cuando se habla de predecir un promedio de \(m\) observaciones asociadas a un valor fijo de \(X_{h}\), entonces la fórmula se simplifica y sólo difiere de la de predicción de un valor por lo que en la varianza se cambia \(1\) por \(1/m\). Lo anterior se puede comprobar viendo que

m=3
var.PredMeanYh_x100v2=MSE*(1/m+1/n+(100-xbar)^2/SSx)
alpha=.1
tcuantil=qt(1-alpha/2, n-2,lower.tail = TRUE)
print(c(beta0+beta1*100-tcuantil*sqrt(var.PredMeanYh_x100v2),
        beta0+beta1*100+tcuantil*sqrt(var.PredMeanYh_x100v2)))
## [1] 365.2356 473.5366

lo cual nos arroja el mismo intervalo de predicción.

\(c)\) Suponga que la compañía ha firmado un contrato para producir en tres de sus unidades, un lote de tamaño 100 en cada uno. El directivo está interesado en saber cuántas horas en total se usarán para cumplir con lo establecido en el contrato.

Solución: Este intervalo se puede obtener multiplicando por tres el obtenido en \(b)\), pues en el inciso anterior trabajamos y hablamos de un promedio, y en este caso se está hablando de las horas totales. Luego

# Este intervalo se puede obtener multiplicando por tres el obtenido en b)

print(3* c(beta0+beta1*100-tcuantil*sqrt(var.PredMeanYh_x100v2),
        beta0+beta1*100+tcuantil*sqrt(var.PredMeanYh_x100v2)) )
## [1] 1095.707 1420.610

o directamente usando la fórmula vista inicilamente. Buscamos pues la estimación de \(\sum_{j=1}^ m y_{hj}\) con \(m=3\) y con \(x_{h1}= x_{h2}=x_{h3}=100\). De donde \(a_j=1\) para \(j=1,2,3\).

sumA_j=(1+1+1)
sumA_j_2=(1^2+1^2+1^2)
SumA_jXhj=((1)*100+(1)*100+(1)*100)
var.PredSumYh_x100=MSE*(sumA_j_2+sumA_j^2/n+(SumA_jXhj-sumA_j*xbar)^2/SSx)
alpha=.1
tcuantil=qt(1-alpha/2, n-2,lower.tail = TRUE)
print(c(sumA_j*beta0+beta1*SumA_jXhj-tcuantil*sqrt(var.PredSumYh_x100),
        sumA_j*beta0+beta1*SumA_jXhj+tcuantil*sqrt(var.PredSumYh_x100)))
## [1] 1095.707 1420.610

En conclusión se usarán entre 1095 y 1420 horas totales para cumplir con lo establecido en el contrato.