Ejemplo práctico

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 estudiar el tamaño óptimo fue a través de 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ó el número de 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=\textrm{tamaño del lote}\)

\(Y=\textrm{horas de trabajo}\)

Los datos son

# Importamos los datos

library(ALSM)

Datos=TolucaCompany

veamos

head(Datos)
##    x   y
## 1 80 399
## 2 30 121
## 3 50 221
## 4 90 376
## 5 70 361
## 6 60 224

Gráficamente podemos ver el diagrama de dispersión

plot(y~x, data=Datos, cex=.9, cex.axis=.7, cex.lab=.8)

Procedamos a efectuar la estimación puntual “a mano”. Recordemos que

\[ \beta_{1}=\frac{\sum_{i=1}^{n}(Y_{i}-\bar{Y})(X_{i}-\bar{X})}{\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}} \]

y que \(\beta_{0}=\bar{Y}-\beta_{1}\bar{X}\). De tal manera tendremos que

options(digits=10)

# cálculos previos

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))

# Estimadores 
beta1=SSxy/SSx
beta0=ybar-beta1*xbar

print(paste("b0:",beta0, "b1:", beta1))
## [1] "b0: 62.3658585858586 b1: 3.57020202020202"

donde la interpretación de \(\beta_{1}\) es que, por cada unidad adicional en un lote, el número promedio de horas requerido aumenta en \(3.57\). De tal manera, la recta ajustada es

\[ \mathbb{E}[y;x]=\hat{\beta_{0}}+\hat{\beta_{1}}x=62.36 + 3.57x \]

Por ejemplo, si queremos estimar el número promedio de horas necesarias cuando se produce un lote de tamaño \(x=65\), bastará con evaluar \(\mathbb{E}[y;65]\):

# Definimos la recta ajustada como una función

model_lineal <- function(x){
  beta0 + beta1 * x 
}

# evaluamos en x=65

model_lineal(65)
## [1] 294.4289899

Por otro lado, los errores observados se obtienen de calcular \(e_{i}=y_{i}-\hat{y_{i}}\). Podemos calcular dicho errores agregando primero una nueva columna al dataFrame de Datos que represente los valores estimados (\(\hat{y_{i}}\))

Datos$yhat = beta0 + beta1 * Datos$x

para despues efectuar la correspondiente diferencia \(y_{i}-\hat{y_{i}}\)

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.577979798

Supongamos que se tienen los siguientes aspectos por responder.

  1. El directivo de la empresa está interesado en saber, con una confianza del 95 %, cuál sería el aumento promedio de horas de trabajo al aumentar en una unidad el tamaño del lote.

  2. ¿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 promedio de trabajo y el tamaño del lote.

  3. 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 %.

  4. 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 3300 dólares por lote?

  5. Suponga que por políticas de la empresa, el número de horas de trabajo por lote se limitará a partir de ahora para que en promedio sea de 400. ¿Cuál sería el tamaño máximo de un lote que se puede producir dada esta limitante?

Solución: i. y ii. En R podemos efectuar de manera rápida el análisis “a mano” que realizamos anteriormente utilizando la función lm().

# realizamos el ajuste 

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

y

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-10

donde los valores de \(\beta_{0}\) y \(\beta_{1}\) se encuentran en coefficientes con Intercept y x respectivamente. Además, este mismo apartado de coefficientes tenemos los valores de \(\sqrt{\hat{V}[\hat{\beta_{0}}]}\) y \(\sqrt{\hat{V}[\hat{\beta_{1}}]}\), los valores de las estadísticas \(t\) asociadas a las pruebas de hipótesis

\[ H_{0}: \beta_{0}=0\ \ vs\ \ H_{a}: \beta_{0}\neq 0 \]

y

\[ H_{0}: \beta_{1}=0\ \ vs\ \ H_{a}: \beta_{1}\neq 0 \]

así como los \(p-values\) correspondientes.

Para calcular intervalos de confianza de los parámetros \(\beta_{0}\) y \(\beta_{1}\) se puede utilizar la función confint()

# intervalo al 95% de confianza 

confint(fit, level=0.95)
##                   2.5 %        97.5 %
## (Intercept) 8.213710748 116.518006423
## x           2.852435427   4.287968613

Con lo anterior ya tenemos la información necesaria para responder las preguntas i. y ii.

Para la i. tenemos que el número promedio de horas aumenta entre 2.8 y 4.28 al aumentar en una unidad el tamaño del lote.

Para la ii. debemos realizar la prueba de hipótesis \(H_{0}: \beta_{1}=0\ \ vs\ \ H_{a}: \beta_{1}\neq 0\). Rechazamos \(H_{0}\) con una significancia de \(\alpha=0.05\) pues \(4.4488e-10<\alpha=0.05\), es decir, podemos decir que hay una relación lineal entre las horas promedio de trabajo y el tamaño del lote.

Para estimar de forma puntual e intervalar \(\mathbb{E}[y;x]\) para ciertos valores de \(x\) se usa la función predict. Con ella podemos realizar, por ejemplo, las prediciones de los datos

# Valores de x
datos_x = Datos$x

# ybarra
datos_predict = as.vector(predict(fit))

# creamos un data frame a partir de los datos en x y ybarra
datos2 = as.data.frame(cbind(datos_x, datos_predict))

datos2
##    datos_x datos_predict
## 1       80   347.9820202
## 2       30   169.4719192
## 3       50   240.8759596
## 4       90   383.6840404
## 5       70   312.2800000
## 6       60   276.5779798
## 7      120   490.7901010
## 8       80   347.9820202
## 9      100   419.3860606
## 10      50   240.8759596
## 11      40   205.1739394
## 12      70   312.2800000
## 13      90   383.6840404
## 14      20   133.7698990
## 15     110   455.0880808
## 16     100   419.3860606
## 17      30   169.4719192
## 18      50   240.8759596
## 19      90   383.6840404
## 20     110   455.0880808
## 21      30   169.4719192
## 22      90   383.6840404
## 23      40   205.1739394
## 24      80   347.9820202
## 25      70   312.2800000

lo cual es justamente lo que se obtuvo al inicio. Podemos calcular un intervalo al 95% de confianza utilizando nuevamente la función predict.

# Definimos un dataframe

valor65 <- data.frame(x=c(65))
valor65
##    x
## 1 65
# Utilizamos la función predict

predValor65 <- predict(fit, valor65, interval="confidence", level=0.95)
head(predValor65)
##           fit         lwr         upr
## 1 294.4289899 273.9129153 314.9450645

donde el número promedio de horas necesarias cuando se produce un lote de tamaño 65 se encuentra entre los valores 273.91 y 314.94, y además \(\mathbb{E}[y,65]=294.4\). Con lo anterior acabamos de responder a la pregunta iii.

Con el paquete multcomp se pueden obtener intervalos y pruebas de hipótesis para las combinaciones lineales de los parámetros, \(\theta=Z_{0}\beta_{0}+Z_{1}\beta_{1}\).

# Cargamos el paquete

library(multcomp)

Por ejemplo, para el inciso III. podemos efectuar, alternativamente, el siguiente código

# Requeriremos una matriz de una fila
matZ0Z1 = matrix(c(1,65), nrow = 1, ncol = 2)

# utilizamos el modelo fit y la especificación de la hipótesis lineal (la matriz matZ0Z1)
prueba2 = glht(fit, linfct = matZ0Z1)

#Intervalo de confianza
confint(prueba2, level = 0.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 294.4289899 273.9129153 314.9450645

lo cual arroja justamente el mismo resultado que en el primer procedimiento.

Ahora bien, tenemos que los ingresos son de 3000 por lote, el costo es el número de horas en promedio trabajadas por el lote 65 multiplicado por 10, esto es \(10\times (\beta_{0}+\beta_{1}65)\). De tal manera, las ganancias son \(3000-10\times (\beta_{0}+\beta_{1}65)\). Podemos considerar la prueba de hipótesis

\[ H_{0}: 3000-10\times (\beta_{0}+\beta_{1}65)\leq 100\ \ vs\ \ H_{a}: 3000-10\times (\beta_{0}+\beta_{1}65)>100 \]

en la cual queda la hipótesis de interés en \(H_{a}\). Equivalentemente tenemos

\[ H_{0}:\beta_{0}+\beta_{1}65\geq 290\ \ vs\ \ H_{a}: \beta_{0}+\beta_{1}65<290 \] luego

# recordemos que líneas arriba definimos la matriz matZ0Z1 

prueba3 = glht(fit, linfct=matZ0Z1, rhs=290, alternative="less")
summary(prueba3)
## 
##   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)

donde la hipótesis de interés se encuentra en la hipótesis alternativa con el símbolo “menor que” o “less that” en inglés, además tendremos two.sided para la igualdad (dada por defecto), greater para \(>\) y less para \(<\).

Notemos de la información arrojada por summary(prueba3) que Pr(<t) es el \(p_value\) asociado a la prueba, de modo que \(p_value=0.67>0.05\). De tal manera, no podemos rechazar \(H_{0}\) con significancia \(\alpha=0.05\), es decir, no hay suficiente evidencia para poder indicar que en promedio se obtendrá una ganancia de 100 por lote cuando se cobran 3000 por cada lote de tamaño 65. Si cobramos 3300 tendremos entonces que contrastar la prueba de hipótesis

\[ H_{0}: \beta_{0}+\beta_{1}65\geq 320\ \ vs\ \ H_{a}: \beta_{0}+\beta_{1}65)<320 \]

y en código

prueba4 = glht(fit, linfct=matZ0Z1, rhs=320, alternative="less")
summary(prueba4)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ x, data = Datos)
## 
## Linear Hypotheses:
##            Estimate Std. Error  t value    Pr(<t)   
## 1 >= 320 294.428990   9.917579 -2.57835 0.0084018 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

donde en este caso \(p_value = 0.0084<0.05=\alpha\). En este caso rechazaremos \(H_{0}\) con significancia de \(\alpha=0.05\), concluyendo en este caso que en promedio se obtendrá una ganancia de más de 100 por lote cuando se cobran 3300 por cada lote de tamaño 65.

  1. Observemos que ahora buscamos \(400=\mathbb{E}[Y; X=x]=\beta_{0}+\beta_{1}x\), de donde

\[ \frac{400-\beta_{0}}{\beta_{1}} \]

es decir, nuestro parámetro de interés es una función no lineal de \(\beta_{0}\) y \(\beta_{1}\). Un estimador puntual se obtendrá usando la propiedad de invarianza de los estimadores máximo verosímilis

\[ \frac{400-\hat{\beta_{0}} }{ \hat{\beta_{1}} } \]

mientras que es posible utilizar la función deltaMethod en la librería car para poder obtener intervalos de confianza asintóticos y aproximados.

library(car)

después

deltaMethod(fit, "(400-b0)/b1", parameterNames = paste("b", 0:1, sep=""), level =0.95)
##                 Estimate         SE      2.5 %    97.5 %
## (400 - b0)/b1 94.5700382  3.6307456 87.4539076 101.68617

el tamaño máximo de los lotes debe ser 95 para cumplir con el límitante de horas.