031 Modelacion para Modelos de Regresion
Eduardo Villarreal
1/5/2021
1 Introducción
Ahora vamos a ver herramientas para modelar distintos casos que puedan presentarse con un ejemplo sobre series de tiempo. Para esto vamos a usar datos de la libreria fpp2 de Rob Hyndman
El objsetivo es predecir la venta del medicamento utilizando regresión.
###Loading Packages -------
#install.packages('fpp2')
require(fpp2)
require(tidyverse)
require(lubridate) #Para manejo de fechasLalibrería fpp2 contiene los datos de una medicina contra la diabetes de nombre a10 que comienza en Julio 1991 y termina en Junio del 2008
Un posible modelo para hacer esta predicción es:
\(y_t = \beta_0 + T_t + S_t+\epsilon_t\)
en donde \(T_t\) es la tendencia en el mes \(t\) y \(S_t\) es la estacionalidad del mes \(t\).
Existen muchas maneras de moelar estos dos componentes. Una forma es:
\(y_t = \beta_0 + \sum_{i=1}^n\beta_i Y_i + \sum_{j=1}^m\phi_jX_j+\epsilon_t\)
en donde \(Y_i\) es el año \(i\) y \(X_j\) es el mes \(j\)
2 Design Matrix
Originalmente los datos pre-cargados a10 son de timo ts; es decir, son un objeto en R de tipo time series. El primer paso es convertir el objeto de series de tiempo en un data frame:
a10_df = data.frame(a10)
head(a10_df)Ahora necesitamos crear la columna de la fecha Date. Sabemos que la serie comienza en Julio 1991 y termina en Junio del 2008 así que podemos crear una secuencia:
a10_df$Date = seq(as.Date("1991-07-1"), as.Date("2008-06-1"), by = "months")
head(a10_df)De esta forma R y en general cualquier lenguaje de programación no entienden que la fecha es una variable que consta de dos componentes: el año y el mes. Entonce, vamos a generar ambas variables con la función year() y month() de lalibreria lubridate
a10_df = a10_df %>%
mutate(year = factor(year(Date)),
month = factor(month(Date)))
head(a10_df)Con esto, vamos a revisar el efecto de ambas variables con un análisis gráfica:
ggplot(a10_df, aes(x = year, y = a10, fill = year)) +
geom_boxplot() +
theme_bw() +
theme(legend.position = 'none') +
labs(title = 'a10 como función del año')## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
La ventaja de tomar en cuenta la variable year como un factor en lugar de una variable numérica es que permite modelar una tendencia ue quizás no es del todo una línea recta.
ggplot(a10_df, aes(x = month, y = a10, fill = month)) +
geom_boxplot() +
theme_bw() +
theme(legend.position = 'none') +
labs(title = 'a10 como función del mes')## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
En la gráfica podemos apreciar que en enero las ventas se incrementan y que en febrero las ventas tienden a disminuir.
La matriz de diseño se ve de la siguiente manera:
x_matrix = model.matrix(a10 ~ year + month, data = a10_df)
head(x_matrix)## (Intercept) year1992 year1993 year1994 year1995 year1996 year1997 year1998
## 1 1 0 0 0 0 0 0 0
## 2 1 0 0 0 0 0 0 0
## 3 1 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0
## 5 1 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 0 0
## year1999 year2000 year2001 year2002 year2003 year2004 year2005 year2006
## 1 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0
## year2007 year2008 month2 month3 month4 month5 month6 month7 month8 month9
## 1 0 0 0 0 0 0 0 1 0 0
## 2 0 0 0 0 0 0 0 0 1 0
## 3 0 0 0 0 0 0 0 0 0 1
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0
## month10 month11 month12
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 1 0 0
## 5 0 1 0
## 6 0 0 1Lo que está pasando es que ne la matriz de diseño se han creado variables artificiales: 1 para cada año y otra para cada mes. Entonces, este modelo tendrá 28 variables y 1 intercepto.
Recordemos que la función lm() ya hace la matriz de diseño en forma automática y no es necesario dar ese input al modelo. Sólo he recreado la matriz de diseño para entender qué esta sucediento en realidad.
3 Estimación del Modelo
Ahora vamos a estimar el modelo:
m1_lm = lm(a10 ~ year + month, data = a10_df)
summary(m1_lm)##
## Call:
## lm(formula = a10 ~ year + month, data = a10_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9668 -0.5498 0.0585 0.6289 3.8130
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5764 0.5427 10.276 < 2e-16 ***
## year1992 1.1001 0.5637 1.952 0.052572 .
## year1993 1.9206 0.5637 3.407 0.000814 ***
## year1994 2.3609 0.5637 4.188 4.45e-05 ***
## year1995 2.9035 0.5637 5.151 6.94e-07 ***
## year1996 3.7427 0.5637 6.640 3.81e-10 ***
## year1997 4.3872 0.5637 7.783 5.92e-13 ***
## year1998 5.0688 0.5637 8.992 3.90e-16 ***
## year1999 6.1228 0.5637 10.862 < 2e-16 ***
## year2000 7.5010 0.5637 13.307 < 2e-16 ***
## year2001 8.6652 0.5637 15.372 < 2e-16 ***
## year2002 9.7151 0.5637 17.235 < 2e-16 ***
## year2003 10.3316 0.5637 18.328 < 2e-16 ***
## year2004 12.6095 0.5637 22.369 < 2e-16 ***
## year2005 13.9309 0.5637 24.714 < 2e-16 ***
## year2006 15.6558 0.5637 27.774 < 2e-16 ***
## year2007 19.4362 0.5637 34.480 < 2e-16 ***
## year2008 20.6657 0.6640 31.121 < 2e-16 ***
## month2 -5.5587 0.3827 -14.524 < 2e-16 ***
## month3 -4.6819 0.3827 -12.233 < 2e-16 ***
## month4 -4.6184 0.3827 -12.067 < 2e-16 ***
## month5 -3.6078 0.3827 -9.427 < 2e-16 ***
## month6 -3.9492 0.3827 -10.319 < 2e-16 ***
## month7 -2.9699 0.3847 -7.720 8.58e-13 ***
## month8 -2.6911 0.3847 -6.995 5.40e-11 ***
## month9 -2.7218 0.3847 -7.075 3.46e-11 ***
## month10 -1.9385 0.3847 -5.039 1.16e-06 ***
## month11 -1.4696 0.3847 -3.820 0.000185 ***
## month12 -0.2243 0.3847 -0.583 0.560707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.116 on 175 degrees of freedom
## Multiple R-squared: 0.9698, Adjusted R-squared: 0.9649
## F-statistic: 200.4 on 28 and 175 DF, p-value: < 2.2e-16De este resultado, podemos obtener que el valor base de las ventas es \(\beta_0 = 5.57\). Si nos encontráramos en el año de 1998 y en el mes de marzo, las ventas estimadas serían \(\hat{y_t} = 5.57 + 5.0688 -4.681 = 5.95\).
Vamos a visualizar los valores predichos vs los valores reales:
par(mfrow = c(1,2))
plot(x = a10_df$a10, y = m1_lm$fitted.values)
plot(x = a10_df$Date, y = m1_lm$fitted.values, type = 'l', col = 'blue')
lines(y = a10_df$a10, x = a10_df$Date, col = 'grey')
4 Predicciones con el modelo
Ahora, supongamos que queremos predecir los siguientes 6 meses con el modelo:
#Crear el data frame para extrapolar hacia adelante nuevas observacines
x_new = data.frame(year = factor(rep(2008, 6)),
month = factor(c(7:12)))
x_new#Hacer la predicción
y_new = predict(m1_lm, x_new)
y_new## 1 2 3 4 5 6
## 23.27221 23.55095 23.52027 24.30354 24.77244 26.01783Si observamos la gráfica de dispersión de los valores reales vs los valores estimados, notamos una curvatura. Esto se debe a que los datos exhiben una tendencia que no es del todo lineal. ¿Cómo corregimos esto?
Recordemos que una forma de “aplanar” los datos (reducir su varianza) es haciendo unatransformación. La más común es tomar el log(y):
\(ln(y_t) = \beta_0 + \sum_{i=1}^n\beta_i Y_i + \sum_{j=1}^m\phi_jX_j+\epsilon_t\)
m1_lm = lm(log(a10) ~ year + month, data = a10_df)
#Graficar
par(mfrow = c(1,2))
plot(x = a10_df$a10, y = exp(m1_lm$fitted.values))
plot(x = a10_df$Date, y = exp(m1_lm$fitted.values), type = 'l', col = 'blue')
lines(y = a10_df$a10, x = a10_df$Date, col = 'grey')
La predicción con este nuevo modelo es:
#Crear el data frame para extrapolar hacia adelante nuevas observacines
x_new = data.frame(year = factor(rep(2008, 6)),
month = factor(c(7:12)))
x_new#Hacer la predicción
y_new = exp(predict(m1_lm, x_new))
y_new## 1 2 3 4 5 6
## 24.30529 24.65124 24.79053 26.57572 27.49932 31.75986Ahora si vemos el gráfico de valores ajustados vs valores reales, la tendencia ya es lineal. Ahora vamos a revisar los supuestos del modelo:
#Análisis de residuales
par(mfrow = c(2,2))
#Normalidad
plot(density(m1_lm$residuals), col = 'red', main = 'Desityplot de residuales') # density plot for 'speed'
polygon(density(m1_lm$residuals), col="red")
abline(v = 0, col = 'blue')
qqnorm(m1_lm$residuals)
qqline(m1_lm$residuals, col = 'red')
#Varianza Constante
plot(m1_lm$residuals, main = 'Constant Variance')
#Residuales independientes
acf(m1_lm$residuals)
En general, los residuale se bien excepto por la gráfica de autocorrelación o ACF en donde podemos ver que algunas de las líneas sobrepasan los límites de confianza, por ejmplo, en los razagos 2, 4, 6, 12 indicando un efecto de autocorrelación de los errores por lo que no se cumple el supuesto de independencia.
Nota: es series de tiempo es muy importante que los residuales sean independientes.Veremos cómo corregir este problema en el tema de series de tiempo. el objetivo de esta lectura es mostrar cómo podemos gregregar variables artificiales a partir de otras que no aparecen de forma exlicita en los datos.
En este link puedes consultar más detalles sobre regresión con series de tiempo https://otexts.com/fpp2/regression.html
5 Usando la librería forecast
La librería forecast fue escrita por Rob Hyndman de la universidad de Monash en Australia y se ha convertido en un estándar para trabajar análisis de series temporales.
El objetivo es mostrar la construcción del modelo de regresión con la librería forecast. Esta librería toma como input un objeto de tipo series de tiempo o ts:
#install.packages('forecast')
require(forecast)La función tslm(formula, data) computa la regresión de seires de tiempo:
#Revisar que a10 es una serie temporal
str(a10)## Time-Series [1:204] from 1992 to 2008: 3.53 3.18 3.25 3.61 3.57 ...#Ajustar el modelo
m1_ts = tslm(log(a10) ~ trend + season, data = a10)
summary(m1_ts)##
## Call:
## tslm(formula = log(a10) ~ trend + season, data = a10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.159457 -0.040126 0.004481 0.040173 0.222543
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.545e+00 1.719e-02 89.877 < 2e-16 ***
## trend 9.369e-03 7.531e-05 124.414 < 2e-16 ***
## season2 -5.212e-01 2.169e-02 -24.030 < 2e-16 ***
## season3 -4.210e-01 2.169e-02 -19.409 < 2e-16 ***
## season4 -4.274e-01 2.169e-02 -19.704 < 2e-16 ***
## season5 -3.480e-01 2.169e-02 -16.042 < 2e-16 ***
## season6 -3.727e-01 2.169e-02 -17.181 < 2e-16 ***
## season7 -2.976e-01 2.169e-02 -13.719 < 2e-16 ***
## season8 -2.928e-01 2.169e-02 -13.500 < 2e-16 ***
## season9 -2.966e-01 2.169e-02 -13.673 < 2e-16 ***
## season10 -2.364e-01 2.169e-02 -10.900 < 2e-16 ***
## season11 -2.116e-01 2.169e-02 -9.757 < 2e-16 ***
## season12 -7.695e-02 2.169e-02 -3.548 0.000489 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06323 on 191 degrees of freedom
## Multiple R-squared: 0.9885, Adjusted R-squared: 0.9878
## F-statistic: 1365 on 12 and 191 DF, p-value: < 2.2e-16La ventaja de esta función esque con el argumento trend y season, automáticamente computa las variables dummies para la tendencia y la estacionalidad sobre el mismo objeto de series de tiempo!. Ya no es necesario computar la fecha para poner estas variables artificiales como hicimos antes.
Podemos visualizar el modelo con autoplot
autoplot(log(a10), series="Data") +
autolayer(fitted(m1_ts), series="Fitted") +
xlab("Year") + ylab("a10") +
ggtitle("Ventas de a10")
Y también podemos obtener un pronóstico que ya computa un intervalo de confianza para la predicción que en pronósticos siempre es útil. Como hemos sacado el logaritmo natural, necesitamos obtener el dato real regresando con exp()
m1_tsfcs = forecast(m1_ts, h = 12)
data.frame(m1_tsfcs) %>% sapply(exp)## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## [1,] 23.75302 21.83217 25.84287 20.87123 27.03270
## [2,] 24.09111 22.14292 26.21071 21.16831 27.41747
## [3,] 24.22724 22.26803 26.35881 21.28792 27.57240
## [4,] 25.97187 23.87158 28.25694 22.82089 29.55792
## [5,] 26.87448 24.70120 29.23897 23.61399 30.58516
## [6,] 31.03821 28.52822 33.76904 27.27257 35.32380
## [7,] 33.83640 31.10013 36.81342 29.73127 38.50834
## [8,] 20.28202 18.64186 22.06649 17.82135 23.08245
## [9,] 22.63078 20.80068 24.62189 19.88515 25.75551
## [10,] 22.69803 20.86249 24.69506 19.94424 25.83205
## [11,] 24.80502 22.79910 26.98744 21.79561 28.22996
## [12,] 24.42727 22.45189 26.57644 21.46368 27.80005#Graficar el pronóstico
autoplot(m1_tsfcs)
#Graficar los residuales del modelo
checkresiduals(m1_ts)
##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals from Linear regression model
## LM test = 62.574, df = 24, p-value = 2.742e-05