Una compañía manufactura equipo de refrigeración, así como partes de reemplazo. En el pasado, una de las partes de reemplazo ha sido producida periódicamente en lotes de varios tamaños. En el momento en el que la compañía realiza un estudio de costos, sus directivos tienen la pregunta de cuál sería el tamaño óptimo de los lotes que se debería producir de esta parte de reemplazo. La producción de esta parte de reemplazo incluye la definición y puesta en marcha del proceso de producción (lo que se debe realizar sin importar el tamaño del lote), el uso de máquinas particulares y operaciones de ensamblado. Una forma de medir el tamaño óptimo fue la relación entre el tamaño del lote y las horas de trabajo total requeridas para su producción. Para esto último se tiene un registro de las últimas 25 producciones en donde se midió las horas de producción y el tamaño del lote. Las condiciones de producción se consideran estables dentro de los 6 meses en las cuales las 25 producciones se llevaron a cabo y se espera que continúen así en los próximos tres años.
Es decir, se tiene lo siguiente
\(X=\) tamaño del lote.
\(Y=\) horas de trabajo.
Los datos son:
#Importando los datos
library(ALSM)
Datos=TolucaCompany
head(Datos)## x y
## 1 80 399
## 2 30 121
## 3 50 221
## 4 90 376
## 5 70 361
## 6 60 224str(Datos)## 'data.frame': 25 obs. of 2 variables:
## $ x: int 80 30 50 90 70 60 120 80 100 50 ...
## $ y: int 399 121 221 376 361 224 546 352 353 157 ...require(ggplot2)
require(ggiraph)
require(ggiraphExtra)
ggPoints(aes(x=x,y=y),smooth=FALSE, data=Datos,interactive=TRUE)¿Qué se observa?
¿Se podría considerar que la relación es lineal?
En el modelo de regresión lineal simple se asume lo siguiente:
\[y_i=\alpha+\beta x_i + \epsilon_i, i=1,..., n,\]
donde \(E(\epsilon_i)=0, \; V(\epsilon_i)=\sigma^2 \;\; \text{y} \;\; Cov(\epsilon_i, \epsilon_j)=0, \; i\neq j \; \; \forall i,j = 1,...,n.\)
En general, se realiza el supuesto adicional siguiente para realizar inferencias sobre los parámetros:
\[y_i \sim N(\mu_i, \sigma^2),\] con \(\mu_i=\alpha+\beta x_i\) y \(Cov(y_i, y_j)=0, \; i\neq j \; \; \forall i,j = 1,...,n.\)
Los estimadores de \(\alpha\) y \(\beta\) obtenidos con ambos métodos coinciden y son:
\[\begin{align} \widehat{\alpha} = \overline{y} - \widehat{\beta}\overline{x} \end{align}\]
\[\begin{align} \widehat{\beta} = \frac{\sum_{i=1}^{n}(x_{i}-\overline{x})(y_{i}-\overline{y})}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} \end{align}\]
Algunas estadísticas que se pueden calcular con los datos son:
\(\sum x_i= 1750\), \(\sum y_i= 7807\), \(\sum y_i^2= 2745173\), \(\sum x_i^2= 142300\), \(\sum x_i y_i= 617180\),
\(\sum (x_i-\overline{x})^2=\sum x_i^2-n\overline{x}^2= 19800\),
\(\sum (y_i-\overline{y})^2=\sum y_i^2-n\overline{y}^2= 307203\), y
\(\sum (x_i-\overline{x})(y_i-\overline{y})=\sum x_iy_i-n\overline{x}\overline{y}= 70690\).
options(digits=10)
(xbar=mean(Datos$x))## [1] 70(SSx=sum((Datos$x-xbar)^2))## [1] 19800(ybar=mean(Datos$y))## [1] 312.28(SSy=sum((Datos$y-ybar)^2))## [1] 307203.04(SSxy=sum((Datos$y-ybar)*(Datos$x-xbar)))## [1] 70690(beta=SSxy/SSx)## [1] 3.57020202(alpha=ybar-beta*xbar)## [1] 62.36585859La interpretación de \(\beta\) es: el número promedio de horas se aumenta por 3.57 horas por cada unidad adicional que se produce en un lote.
La recta ajustada es:
\(\widehat{E(y|x)}= \hat{y}= \hat{\alpha}+\hat{\beta} x= 62.366+3.5702x\)
Por ejemplo, si se quiere estimar el número promedio de horas necesarias cuando se produce un lote de tamaño 65, basta evauar la recta ajusta en x=65.
\(\widehat{E(y|x=65)}= \hat{y}= 62.366+3.5702*65=294.4\)
La recta que se graficó es entonces la recta ajusta evaluada, por ejemplo, para x=80.
\(\widehat{E(y|x=80)}= \hat{y}= 62.366+3.5702*80=347.98\)
Los residuales se obtienen al calcular la diferencia entre el valor observado \(y_i\) y el obtenido al evaluar la recta de regresión \(\hat{y}_i\). Estos los hemos denotado como \(\hat{\epsilon}_i=y_i-\hat{y}_i\).
options(digits=10)
Datos$yhat=alpha+beta*Datos$x
Datos$error=Datos$y-Datos$yhat
head(Datos)## x y yhat error
## 1 80 399 347.9820202 51.017979798
## 2 30 121 169.4719192 -48.471919192
## 3 50 221 240.8759596 -19.875959596
## 4 90 376 383.6840404 -7.684040404
## 5 70 361 312.2800000 48.720000000
## 6 60 224 276.5779798 -52.577979798Un estimador insesgado de \(\sigma^2\) es
\[\begin{align} \widehat{\sigma}^2= \dfrac{\sum_{i=1}^n (y_i-\widehat{\alpha}-\widehat{\beta} x_i)^2}{n-2} = \dfrac{\sum_{i=1}^n (y_i-\overline{y})^2-\widehat{\beta}^2\sum_{i=1}^n (x_i-\overline{x})^2}{n-2}\end{align} \]
n=length(Datos$x)
MSE=sum((Datos$error)^2)/(n-2)
print(MSE)## [1] 2383.715617En general, supongamos que se tiene interés en estimar
\[\theta=z_0 \alpha + z_1 \beta,\]
donde \(z_0\) y \(z_1\) son constantes, no ambas iguales a cero.
Sabemos lo siguiente:
De acuerdo con el Teorema de Gauss-Markov, el estimador \(z_0 \widehat{\alpha} + z_1 \widehat{\beta}\) es el de mínima varianza entre los estimadores lineales e insesgados.
Un estimador insesgado de \(V(\widehat{\theta})\) es \(\hat{V}(\widehat{\theta})= \widehat{Var}(z_0 \widehat{\alpha} + z_1 \widehat{\beta})=\hat{\sigma}^2[\dfrac{z_0^2}{n}+\dfrac{(z_1-z_0\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]\)
Un intervalo de confianza para \(\theta\) es \[[\hat{\theta} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\theta})}, \; \hat{\theta} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \widehat{Var}(\hat{\theta})}] \]
Considerando la siguiente prueba de hipótesis \[H_0: z_0 \alpha + z_1 \beta = c_0 \] vs \[H_a: z_0 \alpha + z_1 \beta \neq c_0 \]
La estadística de prueba es \[t^*= \dfrac{z_0 \widehat{\alpha} + z_1 \widehat{\beta} -c_0}{\sqrt{ \widehat{Var}(z_0 \widehat{\alpha} + z_1 \widehat{\beta})}}\]
La regla de decisión para una significancia \(\delta\) es: Rechazar \(H_0\) si \(|t^*|> t_{n-2, 1-\frac{\delta}{2}}\).
Casos particulares son:
Ejemplos:
El directivo de la empresa quiere saber un intervalo de confianza al 95% para el aumento promedio de horas de trabajo al aumentar en una unidad el tamaño del lote.
¿Se puede decir que el modelo tiene sentido? Es decir, que al usar la variable tamaño del lote se puede decir que hay relación lineal entre las horas de trabajo y el tamaño del lote.
Un directivo está interesado en saber cuántas horas en promedio se requieren para producir un lote de tamaño 65. Dé una estimación intervalar al 95%.
¿Qué significan las bandas en gris de la gráfica anterior? ¿por qué su aplitud es menor alrededor de \(x=70\)?
Suponga que el costo de producción por hora es de 10 dólares. Si se decidierá producir de ahora en adelante lotes de tamaño 65 y cobrar 3000 dólares por lote, los directivos quieren saber si en promedio se puede tener una ganancia de más de 100 dólares por lote. Realizar una prueba de hipótesis. ¿Y si cobrarán 3100 dólares por lote?
Ejemplos.
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. De un intervalo al 90% de confianza. En este caso se usaría un intervalo de predicción para una nueva observación \(y_h\), asumiendo un valor \(x_h\). \[[\hat{\theta} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \hat{\sigma}^2[1+\dfrac{1}{n}+\dfrac{(x_h-\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]}, \; \hat{\theta} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \hat{\sigma}^2[1+\dfrac{1}{n}+\dfrac{(x_h-\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]}] \] Notar que una predicción puntual del valor de \(y_h\) es \(\hat{\theta}=\hat{y}= \hat{\alpha}+\hat{\beta} x_h\).
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 se usarán para cumplir la producción y cumplir con lo establecido en el contrato. En este caso se usaría un intervalo de predicción para el promedio de \(m\) nuevas observaciones \(y_{h1}\), \(y_{h2}\),…, \(y_{hm}\), asumiendo un valor igual de \(x_h\). Es decir, el objetivo es estimar \[\dfrac{\sum_{j=1}^{m}y_{hj}}{m}.\] Un intervalo de predicción para \(\dfrac{\sum_{j=1}^{m}y_{hj}}{m}\) se obtiene como: \[[\hat{y} - t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \hat{\sigma}^2[\dfrac{1}{m}+\dfrac{1}{n}+\dfrac{(x_h-\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]}, \; \hat{y} + t_{n-2, 1-\frac{\delta}{2}}\sqrt{ \hat{\sigma}^2[\dfrac{1}{m}+\dfrac{1}{n}+\dfrac{(x_h-\overline{x})^2}{\sum_{i=1}^{n}(x_{i}-\overline{x})^2} ]}], \] donde \(\hat{y}= \hat{\alpha}+\hat{\beta} x_h\)
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 la producción y cumplir con lo establecido en el contrato.
Usando la tabla anova, ¿qué puede concluir sobre la bondad del modelo propuesto?
¿Qué porcentaje de la variabilidad en las horas de trabajo se explica con el modelo que incluye la variable tamaño de lote?
Sol. Un intervalo para \(\beta\).
var.beta=MSE/SSx
delta=.05
qt(1-delta/2, n-2,lower.tail = TRUE)## [1] 2.06865761print(c(beta-qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.beta),
beta+qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.beta)))## [1] 2.852435427 4.287968613Sol. Prueba de hipótesis \(H_0: \beta=0\) vs \(H_a: \beta \neq 0\). Usando una significancia de \(\delta=.05\).
test=(beta)/sqrt(var.beta)
delta=.05
t= qt(1-delta/2, n-2,lower.tail = TRUE)
print(abs(test))## [1] 10.28959226print(t)## [1] 2.06865761Como \(|test|=10.28959226>2.068658\) entonces se rechaza \(H_0\). Es decir, se puede decir que hay una relación lineal entre las horas de trabajo y el tamaño del lote.
Sol. Calcular un intervalo de confianza para \(E(y|x=65)=\alpha+\beta 65\).
var.Eyx65=MSE*(1/n+(65-xbar)^2/SSx)
delta=.05
qt(1-delta/2, n-2,lower.tail = TRUE)## [1] 2.06865761print(c(alpha+beta*65-qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.Eyx65),
alpha+beta*65+qt(1-delta/2, n-2,lower.tail = TRUE)*sqrt(var.Eyx65)))## [1] 273.9129153 314.9450645Con una confianza de 90%, el número de horas requeridas cuando se producen lotes de tamaño 65 está contenido entre 273.9 y 314.9.
Sol. Visto en clase. Corresponden a los intervalos de confianza para \(E(y|x)\), donde \(x\) son los valores en el eje de las \(x\). Nota. Estos intervalos se leen para cada valor de \(x\) y no de forma simultánea.
Sol. Se debe realiza la prueba de hipótesis \(H_0: \alpha + (65) \beta \ge 290\) vs \(H_a: \alpha + (65) \beta < 290\).
var.Eyx65=MSE*(1/n+(65-xbar)^2/SSx)
test=(alpha+beta*65-290)/sqrt(var.Eyx65)
delta=.05
t= -qt(1-delta, n-2,lower.tail = TRUE)
print(test)## [1] 0.4465797598print(t)## [1] -1.713871528Como \(test\) no es menor a \(t\), entonces no se rechaza \(H_0\), es decir, no hay suficiente evidencia para poder indicar que en promedio se obtendrá una ganancia de más de 100 dólares por lote.
Si se cobrarán 3100 dólares, la prueba de hipótesis \(H_0: \alpha + (65) \beta \ge 300\) vs \(H_a: \alpha + (65) \beta < 300\).
var.Eyx65=MSE*(1/n+(65-xbar)^2/SSx)
test=(alpha+beta*65-300)/sqrt(var.Eyx65)
delta=.05
t= -qt(1-delta, n-2,lower.tail = TRUE)
print(test)## [1] -0.5617308708print(t)## [1] -1.713871528Tampoco se rechaza \(H_0\).
El ajuste en R se realiza usando la función lm.
El summary nos proporciona, los valores de los coeficientes \(\hat{\alpha}\) y \(\hat{\beta}\), así como \(\sqrt(\hat{V}(\hat{\alpha}))\) y \(\sqrt(\hat{V}(\hat{\beta}))\), los valores de las estadísticas t asociados a las pruebas de hipótesis \(H_0: \alpha=0\) vs \(H_0: \alpha\ne 0\) y \(H_0: \beta=0\) vs \(H_0: \beta\ne 0\), y los \(p-values\) correspondientes. Con los p-values, la regla a usar es la siguiente. Si el p-value es menor a la significancia \(\delta\), entonces se rechaza \(H_0\).
Además también nos proporciona \(\hat{\sigma}\), \(n-2\), \(R^2\) y los resultados de la prueba \(F\) de la tabla anova.
fit=lm(y~x, data=Datos)
summary(fit)##
## Call:
## lm(formula = y ~ x, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -83.87596 -34.08808 -5.98202 38.82606 103.52808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.3658586 26.1774339 2.38243 0.025851 *
## x 3.5702020 0.3469722 10.28959 4.4488e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.82331 on 23 degrees of freedom
## Multiple R-squared: 0.8215335, Adjusted R-squared: 0.8137741
## F-statistic: 105.8757 on 1 and 23 DF, p-value: 4.448828e-10Para calcular intervalos de confianza de los parámetros.
confint(fit, level=.95)## 2.5 % 97.5 %
## (Intercept) 8.213710748 116.518006423
## x 2.852435427 4.287968613Para estimar \(E(y|x)\) se usa predict, también nos puede dar intervalos de confianza.
newdata <- data.frame(x = c(65) )
Eydadx65 <- predict(fit, newdata, interval = "confidence", level = 0.95)
head(Eydadx65)## fit lwr upr
## 1 294.4289899 273.9129153 314.9450645Nota: Para predecir el valor de \(y\) dado que se observa \(x_h\) también se usa predict, también nos puede dar intervalos de predicción.
newdata <- data.frame(x = c(65) )
Predx65 <- predict(fit, newdata, interval = "prediction", level = 0.95)
head(Predx65)## fit lwr upr
## 1 294.4289899 191.3676044 397.4903754Con el paquete multcomp se pueden obtener intervalos y pruebas de hipótesis para combinaciones lineales de los parámetros.
library(multcomp)
#Para beta
MatZ0Z1=matrix(c(0,1), ncol=2, nrow=1)
c0=0
prueba1=glht(fit, linfct=MatZ0Z1, rhs=c0)
summary(prueba1)##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = y ~ x, data = Datos)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 3.5702020 0.3469722 10.28959 4.4488e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)### Intervalo de confianza
confint(prueba1, level=.95)##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = y ~ x, data = Datos)
##
## Quantile = 2.0686576
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 3.5702020 2.8524354 4.2879686#Para ejemplo v
MatZ0Z1=matrix(c(1,65), ncol=2, nrow=1)
c0=290
#alternative: "two.sided" (default), "greater" or "less"
prueba2=glht(fit, linfct=MatZ0Z1, rhs=c0, alternative="less")
summary(prueba2)##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = y ~ x, data = Datos)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(<t)
## 1 >= 290 294.428990 9.917579 0.44658 0.67032
## (Adjusted p values reported -- single-step method)