Modelo de Regresión Polinómica no paramétrica

  1. Trabajar con los datos del día elegido en el ejercicio anterior.

  2. Ajustar un modelo de regresión polinómica no paramétrica para los datos del día elegido, usando los niveles de \(NO_2\) como respuesta y la hora del día como variable independiente.

  3. Compare este modelo polinómico con el modelo de series de Fourier que ajustó en la tarea anterior.

Librerías necesarias

library(readxl)
library(dplyr)
library(lubridate)
library(ggplot2)
library(gridExtra)
library(purrr)

Carga de Datos

Para efectos de este ejercicio se hace uso del dataset \(DatosLaFlora\) el cual contiene información sobre mediciones horarias de \(NO_2\) y \(O_3\) en (\(\mu g/m3\)) y sobre algunas variables de clima en el norte de Cali (Colombia). Se dispone además de un campo que informa sobre el día y la hora a la que se tomó la medición.

Datos <- read_excel(ruta)

Ajuste de los datos

Lo que se busca en esta sección es conformar un dataframe con las variables (\(NO_2\), \(dia\) y \(hora\_normada\)) del dataset \(DatosLaFlora\), las variables \(dia\) y \(hora\_normada\) se obtienen a partir de la variable \(Fecha \& Hora\); \(hora\_normada\) se usa para restringirnos al intervalo \([0,1]\), este espacio de las funciones cuadrado integrables se denota en general como \(L_2[0,1]\) el cual satisface las propiedades generales de un espacio vectorial.

La transformación se hace como sigue:

para cada hora \(i \in [1,24]\) se aplica la siguiente ecuación para normalizar:

\[x_i=\frac{(2_i-1)}{2n}, \hspace{1 cm} \forall i=1,2,...,24\]

x <- Datos  %>%
        dplyr::select(`Fecha & Hora`, NO2) %>%
        mutate(dia = day(`Fecha & Hora`)) %>%
        mutate(fecha = as.Date(Datos$`Fecha & Hora`)) %>%
        filter(dia %in% 10) %>%
        
        mutate(hora = row_number()) %>%
        mutate(hora_normada = (2*hora - 1)/(2*24)) %>%
        
        dplyr::select(fecha, hora_normada, NO2)

Obtenemos entonces el siguiente \(tibble\):

## # A tibble: 24 x 3
##    fecha      hora_normada   NO2
##    <date>            <dbl> <dbl>
##  1 2011-08-10       0.0208  63.8
##  2 2011-08-10       0.0625  46.1
##  3 2011-08-10       0.104   47.7
##  4 2011-08-10       0.146   39.9
##  5 2011-08-10       0.188   35.2
##  6 2011-08-10       0.229   31.2
##  7 2011-08-10       0.271   26.8
##  8 2011-08-10       0.312   22.9
##  9 2011-08-10       0.354   37.0
## 10 2011-08-10       0.396   35.3
## # ... with 14 more rows

Exploración de los datos

Generamos un gráfico de puntos con las mediciones de \(NO_2\) vs \(hora\_normada\) para observar el comportamiento de los datos a lo largo del día \(10\).

ggplot()+
  geom_point(data = x, aes(x = hora_normada, y = NO2)) +
  labs(y = expression(NO[2]))

Ajuste del modelo

De acuerdo con el Teorema de Taylor, si tenemos una función continua \(\mu\) en el intervalo \([a, b]\) que tiene \(\lambda+1\) derivadas continuas en el intervalo \((a, b)\), entonces si \(x\) y \(c\) son dos puntos en \((a, b)\), se sigue que \(\mu(x)\) se puede representar como:

\[\mu(x) = \mu(c)+\mu'(c)(x-c)+\frac{\mu''(c)}{2!}(x-c)^2+\dots+\frac{\mu^{(\lambda)}(c)}{\lambda!}(x-c)^\lambda+r_{\lambda}(x)\]

El término \(r_{\lambda}(x)\) se conoce como el \(\lambda-residuo\).

La expresión de \(\mu\) usando el Teorema de Taylor puede reescribirse como:

\[\mu(x)=\sum_{j=1}^\lambda \theta_j \mu_j(x) + r_{\lambda}(x)\]

La función \(\mu\) se puede representar como la suma de un polinomio \(p_\lambda(x)\) de grado \(\lambda\) más un residuo \(r_\lambda (x)\); En la medida que \(r_\lambda (x)\) sea más pequeño, entonces \(p_\lambda(x)\) se parecerá más a \(\mu(x)\), entonces existe un \(\lambda\) tal que:

\[\mu(x)\dot{=}\sum_{j=1}^\lambda \theta_j \mu_j(x)\]

donde \(\theta_1=\mu(c)\), \(\theta_j=\frac{\mu^{j-1}(c)}{(j-1)!}\), \(j=2,\dots,\lambda\) y \(\mu_j(x)=(x-c)^{j-1}\), \(j=1,\dots,\lambda\).

Si asumimos que \(\mu\) es una función de regresión que tiene \(\lambda + 1\) derivadas continuas y puede representarse en la forma presentada, entonces, si \(r_\lambda(x)\) es uniformemente pequeño, podríamos escribir:

\[\mu(x)\dot{=}\sum_{j=1}^\lambda \theta_j \mu_j(x) + \epsilon_j\] Diremos entonces que nuestras observaciones se comportan aproximadamente según un modelo de regresión (no-paramétrica) polinómica, el cual presupone que los residuos \(r_\lambda(x_1),\dots,r_\lambda(x_n)\) de la expansión de \(\mu\) en series de Taylor han sido incorporados en los errores aleatorios del modelo básico.

Para encontrar los coeficientes \(\theta_j\) se podría utilizar el método de mínimos cuadrados.

construimos entonces un \(dataframe\) con nuestra variable respuesta y como predictoras cada evaluación de \(x\), \(\lambda\) y \(c\) dados, en la base polinómica mencionada anteriormente, seleccionamos un \(\lambda=9\) y \(c=0.01\) de forma arbitraria para efectos del ejercicio, por lo tanto tenemos como predictoras \(\mu_1(x)+\mu_2(x)+,\dots,+\mu_9(x)\).

Una vez construido el data frame procedemos con la estimación por mínimos cuadrados de los coeficientes \(\theta_j\) haciendo uso de la función lm().

base.poly <- function(x, lambda, c){
        (x-c)^(lambda-1)
}
#-------------------------------
lm.group <- function(x, lambda, c){
        df <-   dplyr::select(x, NO2)
        hora_normada <- x$hora_normada; i <- list(); p <- list()
        
        for (i in 2:lambda) {
                p[[i]] <- base.poly(hora_normada, i, c)
        }
        
        p <- data.frame(matrix(unlist(p), ncol = lambda-1))
        df <-  bind_cols(df, p)
#-------------------------------
        p.i <- df %>% 
                dplyr::select(contains("x")) %>% 
                colnames()
        
        lambda <- lambda %>%
                tibble(lambda = . )
        
        p.i_sum <- paste(p.i, collapse = "+")
        mdl_formula <- as.formula(paste("NO2", p.i_sum, sep = "~"))
        
        mdl <- lm(mdl_formula, data = df)
#-------------------------------
        fitted <- predict(mdl) %>%
                tibble(fitted = . )
        return(fitted)
}

Obtenemos entonces los fitted values del modelo y los agregamos al final del dataframe x.

lambda <- 9
c <- 0.01
fitted <- lm.group(x, lambda, c)
x <- bind_cols(x, fitted)
## # A tibble: 24 x 4
##    fecha      hora_normada   NO2 fitted
##    <date>            <dbl> <dbl>  <dbl>
##  1 2011-08-10       0.0208  63.8   62.3
##  2 2011-08-10       0.0625  46.1   50.4
##  3 2011-08-10       0.104   47.7   45.2
##  4 2011-08-10       0.146   39.9   39.7
##  5 2011-08-10       0.188   35.2   33.6
##  6 2011-08-10       0.229   31.2   29.0
##  7 2011-08-10       0.271   26.8   27.6
##  8 2011-08-10       0.312   22.9   30.0
##  9 2011-08-10       0.354   37.0   35.2
## 10 2011-08-10       0.396   35.3   41.1
## # ... with 14 more rows

Visualización

Finalmente con el nuevo dataframe x podemos generar un gráfico donde se puedan observar los datos reales (puntos) y los datos ajustados por el modelo (línea) de la variable \(NO_2\) vs \(hora\_normada\).

ggplot()+
        geom_point(data = x, aes(x = hora_normada, y = NO2)) +
        geom_line(data = x, aes(x =hora_normada, y = fitted)) +
        labs(y = expression(NO[2])) +
        labs(title = "Estimación con polinomios de Taylor") +
        labs(subtitle = expression(lambda==9 & C==0.01))

Tenemos entonces la estimación del comportamiento de \(NO_2\) para el dia \(10\) usando una base polinomica no parametrica con un \(\lambda=9\) y un \(c=0.01\) seleccionados arbitrariamente para efectos del ejercicio.

El ajuste de este modelo se llevo acabo basándose en material expuesto en (Olaya, 2012).

Regresion Polinomica no parametrica vs series Fourier

En esta sección buscamos observar gráficamente cuál método captura mejor el comportamiento del \(NO_2\) a lo largo del tiempo, usando el mismo valor de \(\lambda\).

plot1 <- ggplot(data = x, aes(x = hora_normada))+
        geom_point(aes(y = NO2)) +
        geom_line(aes(y = fitted, color = "Base polinomios")) +
        geom_line(aes(y = fitted_tarea_5, color = "Base cosenos")) +
        scale_color_discrete(name='Ajuste') +
        labs(y = expression(NO[2]))

plot2 <- ggplot(data = x, aes(x = hora_normada))+
        geom_point(aes(y = NO2)) +
        geom_line(aes(y = fitted), color = "#00CCFF") +
        labs(y = expression(NO[2])) +
        labs(title = "Estimación con polinomios de Taylor") +
        labs(subtitle = expression(lambda==9 & C==0.01))

plot3 <- ggplot(data = x, aes(x = hora_normada))+
        geom_point(aes(y = NO2)) +
        geom_line(aes(y = fitted_tarea_5), color = "#FF6666") +
        labs(y = expression(NO[2])) +
        labs(title = "Estimación con base de cosenos") +
        labs(subtitle = expression(lambda==9))

grid.arrange(plot2, plot3, ncol=2)

plot1

Comparación del MSE

NO_2 <- pull(x, var = NO2)
MSE.Taylor <- (1/24) * sum((NO_2 - fitted)^2)
MSE.fitted.Fourier <- (1/24) * sum((NO_2 - fitted_tarea_5)^2)

ajuste <- c("Series de Fourier", "Polinomios de Taylor")
MSE <- c(MSE.fitted.Fourier, MSE.Taylor)
comparacion <- data.frame(Ajuste = ajuste, MSE = MSE)

plot4 <- ggplot(data = comparacion, aes(x = reorder(x = Ajuste, X = MSE), 
                               y = MSE, color = Ajuste, 
                               label = round(MSE,2))) + 
        geom_point(size = 10) + 
        geom_text(color = "white", size = 2) + 
        labs(x = "Modelo regresión", y = "MSE") + theme_bw() + 
        coord_flip() + theme(legend.position = "none") +
        labs(title = "Comparacion de Modelos") +
        theme(plot.title = element_text(hjust = 0.5))
plot4

Dado el MSE de ambos modelos podriamos decir que la estimación con series de Fourier captura mejor el comportamiento de los datos en este caso, pero teniendo en cuenta que la estimación con polinomios de Taylor no se realizo con ningun criterio de elección del mejor \(\lambda\) no podríamos concluir sobre cual es el mejor método de estimación para este caso, puesto que no sabemos si \(\lambda=9\) y \(c=0.01\) son las mejores elecciones para el metodo de estimación con polinomios de Taylor.

Bibliografía

  • Olaya, J. (2012). Métodos de Regresión No Paramétrica. Universidad del Valle.

  • R Core Team. (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. http://www.r-project.org/