En la práctica de la gestión (aunque no solo en la industria) a menudo es necesario poder anticipar cuál será el resultado de una determinada situación para poder actuar en previsión.
Una vez decidido un modelo, el ajuste por el método de mínimos cuadrados ordinarios se realiza usando la notación de fórmulas y la función lm.
Una fórmula en R indica una relación lineal entre variables o combinaciones de variables.
y ~ xy ~ x - 1 o y ~ 0 + xy ~ x1 + x2y ~ (X1 + x2) ^ 2y ~ x + I(x ^ 2) o y ~ poly(x, 3)y ~ exp(x) o log(y) ~ log(x)lmx <- seq(1, 2, length.out = 11)
y <- round(3 * x + 5 + rnorm(11) / 10, 3)
datos <- data.frame(x = x, y = y)| x | y |
|---|---|
| 1.0 | 7.944 |
| 1.1 | 8.277 |
| 1.2 | 8.756 |
| 1.3 | 8.907 |
| 1.4 | 9.213 |
| 1.5 | 9.672 |
| 1.6 | 9.846 |
| 1.7 | 9.973 |
| 1.8 | 10.331 |
| 1.9 | 10.655 |
| 2.0 | 11.122 |
##
## Call:
## lm(formula = y ~ x, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14286 -0.06881 -0.01278 0.06913 0.15418
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0325 0.1495 33.67 8.87e-11 ***
## x 2.9902 0.0975 30.67 2.04e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1023 on 9 degrees of freedom
## Multiple R-squared: 0.9905, Adjusted R-squared: 0.9895
## F-statistic: 940.5 on 1 and 9 DF, p-value: 2.041e-10
| x | y | yc |
|---|---|---|
| 1.0 | 7.944 | 8.022727 |
| 1.1 | 8.277 | 8.321746 |
| 1.2 | 8.756 | 8.620764 |
| 1.3 | 8.907 | 8.919782 |
| 1.4 | 9.213 | 9.218800 |
| 1.5 | 9.672 | 9.517818 |
| 1.6 | 9.846 | 9.816836 |
| 1.7 | 9.973 | 10.115854 |
| 1.8 | 10.331 | 10.414873 |
| 1.9 | 10.655 | 10.713891 |
| 2.0 | 11.122 | 11.012909 |
Un ajuste es adecuado si
##
## Call:
## lm(formula = y ~ x, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14286 -0.06881 -0.01278 0.06913 0.15418
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0325 0.1495 33.67 8.87e-11 ***
## x 2.9902 0.0975 30.67 2.04e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1023 on 9 degrees of freedom
## Multiple R-squared: 0.9905, Adjusted R-squared: 0.9895
## F-statistic: 940.5 on 1 and 9 DF, p-value: 2.041e-10
##
## Shapiro-Wilk normality test
##
## data: ajuste$residuals
## W = 0.93572, p-value = 0.4716
Si el modelo es correcto (y solo en este caso), este puede usarse para estimar valores de la variable dependiente.
## 1 2
## 8.381549 8.710469
A esta predicción puede añadirse intervalos de confianza (por defecto, a un nivel de confianza del 95 %). Estos se calculan de forma distinta cuando se quiere predecir el valor medio (mejor estimación) o un nuevo dato.
Si x es 1,12, la estimación del promedio es
## fit lwr upr
## 1 8.381549 8.272507 8.490592
y la estimación para un nuevo dato es
## fit lwr upr
## 1 8.381549 8.125803 8.637295
A partir de los datos que figuran en el archivo “cu.txt”, determina el coeficiente de expansión térmica del cobre (k) a las temperaturas (T) de 600 K y de 200 K.
http://www.itl.nist.gov/div898/handbook/pmd/section6/pmd641.htm
Se denomina serie temporal a una secuencia de puntos indexados temporalmente. Es decir, obtenidos a intervalos de tiempo constantes.
En estos casos, el índice temporal (o posición) se incorpora como variable de modelización. El tiempo no se usa como variable independiente.
Solo se presentará el tratamiento del ajuste de serie temporales completas.
Sean las ventas de un determinado producto entre enero de 2013 y mediados de julio de 2016 las que figuran en el archivo “ventas.csv”, ¿cuál seria una previsión razonable de las ventas de julio y agosto de 2016?
| ano | mes | ventas |
|---|---|---|
| 2013 | 1 | 3038 |
| 2013 | 2 | 3420 |
| 2013 | 3 | 3138 |
| 2013 | 4 | 4137 |
| 2013 | 5 | 4302 |
| 2013 | 6 | 3374 |
y ~ 1 )y ~ x )y ~ 0 + x )También pueden usarse funciones no lineales (polinomios, logaritmos, exponenciales…).
## fecha fit lwr upr
## 43 2016-07 3331.548 1740.936 4922.159
## 44 2016-08 3331.548 1740.936 4922.159
## fecha fit lwr upr
## 43 2016-07 2976.685 1364.645 4588.726
## 44 2016-08 2960.180 1342.917 4577.443
La estacionalidad se incorpora a los modelos a través de variables ficticias
datosAj <- datos[1:42,]
datosAj$i <- 1:42
datosAj <- rbind(datosAj,
list(ano=c(2016,2016),
mes=7:8,
ventas=c(NA,NA),
fecha=c("2016-07","2016-08"),
i=43:44))
datosAj$mes <- factor(datosAj$mes)
ajuste <- lm(ventas~i+mes, datosAj[1:42,])
datosAj <- cbind(datosAj,
predict(ajuste, datosAj, interval="prediction"))## fecha fit lwr upr
## 43 2016-07 3996.2143 3373.7414 4618.687
## 44 2016-08 932.2143 309.7414 1554.687
##
## Call:
## lm(formula = ventas ~ i + mes, data = datosAj[1:42, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -463.64 -101.52 -22.13 128.66 538.91
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3071.358 141.707 21.674 < 2e-16 ***
## i -18.769 3.274 -5.732 3.33e-06 ***
## mes2 404.769 180.088 2.248 0.032375 *
## mes3 530.038 180.177 2.942 0.006356 **
## mes4 737.557 180.326 4.090 0.000313 ***
## mes5 785.575 180.534 4.351 0.000153 ***
## mes6 736.344 180.801 4.073 0.000328 ***
## mes7 1731.917 194.485 8.905 8.56e-10 ***
## mes8 -1313.314 194.512 -6.752 2.07e-07 ***
## mes9 1593.788 194.595 8.190 4.96e-09 ***
## mes10 1171.557 194.733 6.016 1.52e-06 ***
## mes11 750.325 194.925 3.849 0.000601 ***
## mes12 1098.761 195.173 5.630 4.42e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 254.6 on 29 degrees of freedom
## Multiple R-squared: 0.9243, Adjusted R-squared: 0.893
## F-statistic: 29.51 on 12 and 29 DF, p-value: 5.199e-13
##
## Shapiro-Wilk normality test
##
## data: ajuste$residuals
## W = 0.98166, p-value = 0.7257
Presentan autocorrelación, es decir dependencia de los valores anteriores.
Consiste en pronosticar el siguiente valor de la serie temporal como igual al último valor disponible en la misma
\[ \\ \hat y_{n+1} = y_n \]
Consiste en pronosticar el siguiente valor de la serie temporal como igual al último valor disponible en la misma más un término dependiente de la variación entre los dos valores anteriores.
\[ \\ \hat y_{n+1} = y_n + k (y_n - y_{n-1}) \]
Si k vale 1, consiste en el ajuste de una recta a los dos datos anteriores.
Consiste en pronosticar el siguiente valor de la serie temporal como media de los m valores anteriores.
\[ \\ \hat{y}_{n+1} = \sum\limits_{i=n-m+1}^{n} \frac{y_i}{m} \]
El ajuste se puede comprobar aplicando el método a los valores anteriores de la serie y calculando los residuales correspondientes.
El análisis de residuales sigue siendo válido en estos casos.
Si los residuales estan normalmente distribuidos, el intervalo de predicción puede calcularse a partir de la desviación estándar de los residuales.
\[ \\ y = \hat y \pm z_{1-\alpha/2} \ s_{res} \]
Disponemos de los precios históricos (promedios mensuales) de uno de los productos de partida de nuestro proceso de fabricación, ¿podemos estimar el precio que nos costará este recurso en el siguiente mes?
Los datos figuran en el archivo “coconutOil.txt”.
Para establecer la calidad de un modelo y tener criterios que permitan comparar modelos, se usan distintos indicadores