Tema 06. Ejercicio 1

Autor/a

Pedro Albarrán

Datos

El siguiente conjunto de datos tiene información sobre el volumen de usuarios de un camino ciclista (“via verde”) en EE.UU. Tenéis más información en la ayuda de RStudio.

library(tidyverse)
library(mosaicData)
data("RailTrail")

La comisión gestora, PVPC, quiere entender la relación entre el volumen de usuarios y variables explicativas como incluyendo la temperatura, lluvia, nubosidad y el día de la semana. Para esto, estimamos el siguiente modelo de regresión:

\[ volume = \beta_0 + \beta_1 \cdot hightemp + \beta_2 \cdot cloudcover + \beta_3 \cdot weekday + \beta_4 \cdot precip + \varepsilon \tag{1}\]

modelo1 <- lm(data = RailTrail,
              volume ~ hightemp + cloudcover + weekday + precip)

Apartado 1

Para este apartado, trabajaréis con un coeficiente concreto de la Ecuación 1 estimada anteriormente. El coeficiente dependerá de la última cifra de vuestro DNI o similar:

  • si es 1, 4 o 7, con el de hightemp
  • si es 2, 5 o 8, con el de cloudcover
  • si es 3, 6 o 9, con el de weekday
  • si es 0, con el de precip
  1. Suponiendo que los errores del modelo siguen una distribución normal, \(\widehat{\beta} \sim N\left(\beta, Var(\widehat{\beta})\right)\) y \((n-k)*s^2 / Var(\varepsilon) \sim \chi^2_{(n-k)}\), donde \(n\) es el número de observaciones, \(k\) es el número de coeficientes (incluida la constante) y \(s^2\) es la estimación de la varianza del error. Calcular el intervalo de confianza al 95% para vuestro coeficiente y el intervalo de confianza al 95% para la varianza del error.
  1. Usar bootstrap para obtener el intervalo de confianza al 95% para vuestro coeficiente y el intervalo de confianza al 95% para la varianza del error. Debéis fijar como semilla vuestro DNI para realizar un bucle como el visto en clase.
  1. Comentar BREVEMENTE las diferencias en los intervalos de confianza de ambos apartados.

Notas

  • Los resultados de una estimación están almacenados en el objeto creado aplicando summary() a la función lm().
sum.modelo1 <- summary(modelo1)

## Para todos los coeficientes (filas) 
## el valor estimado y su error estándar (dos columnas)
sum.modelo1$coefficients[,1:2]

## Varianza del error del modelo
sum.modelo1$sigma**2 

## R-cuadrado
sum.modelo1$r.squared
  • Intervalo de confianza al 95% para el coeficiente estimado: \(\widehat{\beta} \pm z_{0.975} \cdot \sqrt{\operatorname{Var}(\widehat{\beta})}\), donde \(z_{0.975}\) es el valor crítico de la distribución normal estándar .

  • Intervalo de confianza al 95% para la varianza del error: \(\left( \frac{(n - k) \cdot s^2}{\chi^2_{0.975, (n - k)}}, \frac{(n - k) \cdot s^2}{\chi^2_{0.025, (n - k)}} \right)\), donde \(\chi^2_{0.975, (n - k)}\) y \(\chi^2_{0.025, (n - k)}\) son los valores críticos de la distribución \(\chi^2\) con \(n - k\) grados de libertad, correspondientes a los percentiles del 97.5% y 2.5%, respectivamente.

  • En R, los valores críticos se obtienen con qnorm() y qchisq() (mirad la ayuda).

Apartado 2

Realizamos un análsis exploratorio de los datos y encontramos la siguiente forma para la relación no lineal entre el número de visitantes y la temperatura.

Presentar en una tabla los resultados de estimar la Ecuación 1 y añadir primero la temperatura al cuadrado, después también la temperatura al cubo, la temperatura elevada a la cuarta potencia y finalmente elevada a la quinta potencia. Comentar qué modelo elegirías, es decir, qué grado del polinomio en temperatura captura mejor la relación no lineal descrita anteriormente. ¿Podríamos tener una relación no lineal distinta de la descrita por un polinomio?

Notas sobre tablas de resultados en Quarto

  • Podemos combinar tidy() (de la biblioteca broom) con las funciones de la biblioteca kableExtra para incluir tablas de estadísticos descriptivos o de resultados de regresión. En el documento de Quarto, se debe incluir la opción results: markup.
library(broom)
library(kableExtra)
modelo1 %>% tidy() %>% kbl() %>% kable_classic()
  • Notad que tras usar tidy() tenemos un conjunto de datos. Por tanto, podemos usar comandos conocidos para manipular la tabla, p.e., no mostrar todas las columnas.
modelo1 %>% tidy() %>% select(term:std.error) %>% kbl() %>% kable_paper()
  • Con modelsummary(), podemos mostrar los resultados de uno o varios modelos en la misma tabla. También debemos incluir la opción results: markup.

Consideramos el siguiente modelo: \[ volume = \beta_0 + \beta_1 \cdot hightemp + \beta_2 \cdot cloudcover + \beta_3 \cdot weekday + \beta_4 \cdot precip + \beta_5 \cdot spring + \beta_6 \cdot summer + \varepsilon \]

library(modelsummary)
modelo2 <- lm(data = RailTrail,
                volume ~ hightemp + cloudcover + weekday + precip + 
                          spring + summer)

modelsummary(list("Modelo 1" = modelo1, "Modelo 2" = modelo2), 
             gof_map = c("nobs", "r.squared", "adj.r.squared", "F", "rmse") ,
              stars = T)
Modelo 1 Modelo 2
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 76.826 45.483
(62.710) (79.724)
hightemp 5.549*** 5.696***
(0.781) (1.164)
cloudcover -8.503* -8.353*
(3.339) (3.435)
weekdayTRUE -35.045 -37.062+
(21.752) (22.280)
precip -106.483* -98.904*
(40.720) (42.137)
spring 31.239
(32.082)
summer 9.424
(47.504)
Num.Obs. 90 90
R2 0.500 0.509
R2 Adj. 0.476 0.473
F 21.208 14.322
RMSE 89.67 88.85

Apartado 3

Realizamos un nuevo análisis exploratorio para la relación entre el número de visitas y la nubosidad (como porcentaje de cielo cubierto por nubes, en una escala continua de 0 a 10).

Dado lo que observamos en el gráfico, vamos a discretizar cloudcover. Primero, consideramos solo dos grupos: hasta 7.5 y más de 7.5. Luego, consideramos tres rangos: entre 0 y 5, entre 5 y 7.5, y mayor de 7.5. Finalmente, consideramos cuatro categorías: [0,2.5], (2.5,5], (5, 7.5] y (7.5, 10]. Presentar en una tabla el modelo de la Ecuación 1, con la variable cloudclover, y todas las variantes donde la hemos discretizado. Discutir qué especificación preferís y por qué.

Notas

  • Como hemos visto en las transparencias, se puede discretizar una variable con cut() (o cut_width(), cut_interval(), etc.), generando un factor con categorías dadas por los puntos de corte: ej., cut(cloudcover, breaks=c(0, 7.5, 10), include.lowest = T)

  • También se podría generar una variable binaria para cada categoria con ifelse() (o if_else())

Apartado 4

Finalmente, vamos a considerar otra variante de la ecuación Ecuación 1 donde los efectos de la temperatura (hightemp) y de la nubosidad (cloudcover) no son constantes, sino que su efecto es heterogéneo en función de otros factores, en este caso si es un día laborable (weekday). Estimaremos el siguiente modelo

\[ \begin{aligned} volume &= \beta_0 + \beta_1 \cdot hightemp + \beta_2 \cdot cloudcover + \beta_3 \cdot weekday + \beta_4 \cdot precip \\ &+ \beta_5 \cdot hightemp \cdot weekday + \beta_4 \cdot cloudcover \cdot weekday + \varepsilon \end{aligned} \tag{2}\]

modelo1.H <- lm(data = RailTrail, volume ~ (hightemp + cloudcover)*weekday + precip)

Presentar en una tabla los resultado de estimar la ecuación Ecuación 1 y de la ecuación Ecuación 2. Discutir si evidencia de que la temperatura y la nubosidad afectan de manera diferente a las visitas en función de otros factores. ¿Qué modelo preferiría?

Entrega

Rellenad este FORMULARIO con vuestros datos y subid

  • vuestro archivo de .qmd

  • el resultado de renderizarlo: bien un archivo autocontenido .html (o .pdf o .docx) o bien un archivo .html y el directorio relacionado con el mismo nombre; en ambos casos, se recomienda comprimir todo para enviarlo.

IMPORTANTE: el nombre de los ficheros que subáis DEBE seguir el siguiente formato que incluye vuestro número de DNI: ej.,

  • Tema06ej1_123456789.qmd

  • Tema06ej1_123456789.zip