Curso de Análisis estadístico en R

La siguiente entrada contiene un curso de Métodos Estadístico de Investigación en R, usando los fundamentos del tidyverso para la manipulación y comunicación de datos. 

Parte I. Análisis exploratorio en R

Clase 1.0. pdf. video. Introducción a R

Clase 1.1. pdf. video. Estructura de datos (tomado del curso de Derek Corcoran)

Clase 2.1. pdf. video. Análisis exploratorio: posición central y dispersión

Clase 2.2. pdf. video. Análisis exploratorio: posición (quintiles), distribuciones y asociación

Clase 2.3. pdf. video. Análisis exploratorio: datos perdidos y valores atípicos

Clase 2.4. pdf. video. Tablas en R

Parte II. Manipulación y comunicación de datos en R

Clase 3.1. pdf. video. Manipulando datos usando dplyr I (select, filter, mutate)

Clase 3.2. pdf. video. Manipulando datos usando dplyr I (summarize, by_group, accros, if, at)

Clase 3.3. pdf. video. Ejercicios usando dplyr

Clase 4.1. pdf. video. Comunicación de datos con ggplot2 y rmakdowm

Clase 4.2. pdf. video. Ejercicios con ggplot2

Clase 5.1. pdf. video. Gestión de datos relacionales en R (join, reshape, long, wide)

Parte III. Inferencia estadística y nociones de ciencia de datos en R

Clase 6.1. pdf. video. Repaso de inferencia estadística

Clase 6.2. pdf. video. Análisis de independencia estadística en R

Clase 6.3. pdf. video. Ejercicios aplicados a independencia y datos relacionales 

Clase 7.1. pdf. video. Métodos de clasificación en R (Statistical Learning)

Clase 8.1. pdf. video. Modelos de regresión y aprendizaje supervisado

Parte IV. Análisis de series temporales

Clase 9.1.   pdf. video. Exploración y agregación

Clase 10.1. pdf. video. Estacionariedad y modelos ARIMA

Clase 10.2. pdf. video. Ejericicios

Clase 11.1. pdf. video. Pronóstico

Clase 12.1. pdf. video. Modelos de volatilidad (GARCH)

Clase 12.2. pdf. video. Markov GARCH

Test de independencia estadística en R

 1. Análisis preliminar

En la siguiente entrada se muestran algunos ejemplos de cómo realizar test de independencia estadística en R, para responder a uno de los ejercicios de la asignatura de Seminario de Investigación (clase 6 de la parte de R). Los ejercicios se realizan con la base utown del libro Principles of Econometrics 4e, by Carter Hill, William Griffiths, and Guay Lim. Los mismos se encuentran disponibles en el paquete PoEdata.

library(PoEdata)

library(tidyverse)
             
data(utown)        
head(utown,5)        

##     price  sqft age utown pool fplace
## 1 205.452 23.46   6     0    0      1
## 2 185.328 20.03   5     0    0      1
## 3 248.422 27.77   6     0    0      0
## 4 154.690 20.17   1     0    0      0
## 5 221.801 26.45   0     0    0      1

1. Análisis de independencia

1.1. Usando ggplo2, presente un histograma de precios y muestre un summary (min, max,median, q1, q5) condicionado a si la casa tiene o no piscina (pool).

 theme_set(theme_minimal())

#histograma
utown%>%
ggplot(aes(price)) +
geom_histogram(col = "black",fill = "darkblue", alpha = 0.6)

# Summary de precios por grupo
tapply(utown$price, utown$pool, summary)

## $`0`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   134.3   214.8   244.3   246.5   276.8   345.2
##
## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   160.3   219.5   249.1   252.0   285.4   340.7

1.2. Obtenga un boxplot donde relacione el tener piscina o no, con el precio de la casa, y establezca la independencia entre las variables.

utown%>%
  ggplot(aes(y=price, x=factor(pool), fill=factor(pool))) +
  geom_boxplot(col = "black", alpha = 0.6)

1.3. A partir de la comparación de promedios, que conclusión podemos obtener y que limitaciones puede enfrentar esta conclusión sobre la posible obtención de un efecto causal.

Al considerar los intervalos del gráfico anterior se puede observar que los promedios son relativamente parecidos, sin embargo, esta comparación no se puede interpretar como un efecto causal, dado que esta simple comparación omite otros factores relevantes que inciden sobre el precio de las casas.

1.4. Presente una tabla de frecuencia, donde los porcentajes se presenten en relación a las tablas, interprete sus resultados y posteriormente realice el test chi2 para determinar la posible independencia entre las variables.

Necesitamos una tabla de frecuencia de precio con pool, pero price es una variable continua, por lo que, debemos discretizarla, antes de realizar la tabla de frecuencia. Una alternativa es crear rangos de precios (cut(price, breaks = 5)), utilizar los quintiles (ejemplo debajo), o separarla en grupos alrededor de alguna medida como la media (price>mean(price)).

# discretizar la variable precio
utown <- utown %>%
  mutate(quintiprecio = ntile(price, 5))

utown %>%
select(quintiprecio, pool) %>%
table() %>%
prop.table(.,2) 

##             pool
## quintiprecio         0         1
##            1 0.2085427 0.1666667
##            2 0.2010050 0.1960784
##            3 0.1959799 0.2156863
##            4 0.2022613 0.1911765
##            5 0.1922111 0.2303922

Los valores son bastante diferentes, eso es evidencia de independencia, pero se debe realizar el test chi2 de independencia entre variables discretas (p.valor=0.5508, no se rechaza ho de independencia):

# discretizar la variable precio
utown %>%
select(quintiprecio, pool) %>%
table() %>%
chisq.test()

##
##  Pearson's Chi-squared test
##
## data:  .
## X-squared = 3.0422, df = 4, p-value = 0.5508

2. Análisis de independencia: Correlación

2.1. Sabiendo que tenemos una medida del size del terreno de la casa (sqft), crear los quintiles de esta variable y muestre una tabla resumen donde indique el precio max, min, mediano, q1 y q3 del precio de las casas para cada quintil.

utown %>%
  mutate(quintiprecio = ntile(price, 5)) %>%
  select(price,quintiprecio) %>%
  group_by(quintiprecio) %>%
  summarise(min = min(price),
    q1 = quantile(price,0.2),
    mean = mean(price),
    q3 = quantile(price,0.6),
    max = max(price))

## # A tibble: 5 x 6
##   quintiprecio   min    q1  mean    q3   max
##          <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1            1  134.  178.  190.  197.  208.
## 2            2  209.  214.  222.  225.  235.
## 3            3  235.  239.  246.  249.  258.
## 4            4  258.  263.  272.  274.  287.
## 5            5  287.  293.  309.  311.  345.

2.2. Obtenga un scatter donde relacione el size de la casa con el precio, interprete los resultados, coloreando los puntos según la casa tenga o no piscina (pool).

utown %>%
ggplot(aes(sqft, price, color =factor(pool)))+
geom_point()+
geom_smooth()

2.3. Realice el test de correlación de pearson y determine si rechaza o no ho.

Se rechaza ho.

cor.test(utown$price,utown$sqft, method = "pearson")

##
##  Pearson's product-moment correlation
##
## data:  utown$price and utown$sqft
## t = 23.367, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5530750 0.6333235
## sample estimates:
##       cor
## 0.5946785

2.4. Utilice el qqplot para realizar un test de normalidad sobre el price y sqft.

car::qqPlot(utown$price)

2.5. En base al test de normalidad, que podemos concluir sobre las medidas de resumen y el análisis de correlación realizado.

Como para las colas se rechaza que la variable se comporte como una normal, se entiende que es preferible usar alguna estimación no parámetrica para obtener el coeficiente de correlación, tal como el método de Spearman.

3. Test de media

3.1. Realice un test de media para testear si las casas con chimeneas son estadísticamente más grandes que aquellas casas que no tienen.

t.test(sqft~fplace,
       data=utown,
       var.equal=F,
       alternative = "greater",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  sqft by fplace
## t = -3.1236, df = 989.19, p-value = 0.9991
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.8778771        Inf
## sample estimates:
## mean in group 0 mean in group 1
##        24.91187        25.48674

3.2. Realice boxplot donde compare el precio de la casa según el size esté por encima o por debajo de la mediana.

utown%>%
  mutate(sqftmedian = sqft>median(sqft)) %>%
  ggplot(aes(y=price, x=factor(sqftmedian), fill=factor(sqftmedian))) +
  geom_boxplot(col = "black", alpha = 0.6)

3.3. Realice un test de media para testear si las casas con edad por encima de la mediana, tienen a tener un precio menor al resto de casa.

utown <- utown %>%
  mutate(medianEdad = as.numeric(age>median(age)))

t.test(price~medianEdad,
       data=utown,
       var.equal=F,
       alternative = "less",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  price by medianEdad
## t = 1.2334, df = 994.12, p-value = 0.8911
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf 7.685081
## sample estimates:
## mean in group 0 mean in group 1
##        249.2916        246.0001

Trabajar datos diarios agrupados por semanas en R: collapse y pronósticos

 En el siguiente ejemplo se utilizan datos diarios de ventas, simulados para presentar algunos ejemplos de con trabajar con series temporales utilizando el tidyverso [datos]. Puntualmente, se muestra como convertir una data en frecuencia diaria a semanal, y como realizar pronósticos con modelos arima usando la data en frecuencia diaria y en frecuencia semanal (vea dplyr como prerrequisito).

 library(readxl)

library(tidyverse)

library(tseries)

library(forecast)

data_simulada_venta <- read_excel("data_simulada_venta.xlsx")

data_simulada_venta

# A tibble: 1,356 x 2

   fecha               ventas

   <dttm>               <dbl>

 1 2016-01-01 00:00:00  0.343

 2 2016-01-02 00:00:00  1.28 

 3 2016-01-04 00:00:00  1.59 

 4 2016-01-05 00:00:00  1.83 

 5 2016-01-06 00:00:00  2.29 

Convertir la data diaria en datos semanales

El collapse de la data es una operación usada para agregar data, por ejemplo, pasar de una data a nivel de individuos a datos a nivel de viviendas, o en el caso se series temporales pasar de datos de unafrecuencia mayor a una menor, por poner un caso, de mensual a trimestral. En este sentido, una de las operaciones más frecuente es agrupar los datos diarios en semanas, para esto hacemos un by_group por cada semana de los distintos años disponibles en nuestros datos luego de generar dos nuevas variables para identificar los años y las semanas.

a.       Creamos las variables (en caso de no estar disponibles) sobre las que deseamos agregar la data. En nuestros casos, usamos el año (anio) y el día de la semana, porque queremos agregar por semana para cada año.

b.       Usamos by_group para agregar la data, combinado con el summarize, que es donde se indica cual es el proceso por el cual se va agregar la data. En este ejemplo, corresponde a la suma de las ventas diarias. Note que aquí pudiésemos utilizar otras funciones, tales como la media, la varianza o el primer día de la semana, para obtener otra agregación de la data:

data_semanal <- data_simulada_venta  %>% 

                mutate(anio = format(fecha, "%Y"),

                       dia_semana_num = format(fecha, "%W")) %>% 

                        group_by(anio,dia_semana_num) %>% 

                        summarise(suma_venta = sum(ventas),

                                       fecha = max(fecha))

data_semanal

# A tibble: 229 x 4

# Groups:   anio [5]

   anio  dia_semana_num suma_venta fecha              

   <chr> <chr>               <dbl> <dttm>             

 1 2016  00                   1.62 2016-01-02 00:00:00

 2 2016  01                  16.6  2016-01-09 00:00:00

 3 2016  02                  40.9  2016-01-16 00:00:00

 4 2016  03                  55.8  2016-01-23 00:00:00

 5 2016  04                  68.3  2016-01-30 00:00:00

 6 2016  05                  88.9  2016-02-06 00:00:00

 7 2016  06                 110.   2016-02-13 00:00:00

 8 2016  07                 120.   2016-02-20 00:00:00

 9 2016  08                 132.   2016-02-27 00:00:00

10 2016  09                 145.   2016-03-05 00:00:00

# ... with 219 more rows

Obtener una data semanal de un solo día de la semana

Otra forma de agregar datos, es obtener una observación (correspondiente a un día en nuestro ejemplo) como representativo de la data agregada (semanal). Aquí, usamos la función wday del paquete lubridate, que nos permite identificar el día de la semana, para posteriormente, obtener un collapse de la data diaria, pero en vez de agregar la data para obtener las semanas (sumamos todos los días para obtener la data semanal), seleccionados uno de los datos semanales. Puntualmente el viernes, pero en caso de otros días de las semanas, solo es necesario modificar el día de la semana seleccionado en el filter.

library(lubridate)

# todos los viernes del primer mes de 2016 o 2017

data_simulada_venta  %>% 

  mutate(anio = format(fecha, "%Y"),

         mes = format(fecha, "%B"),

         dia_semana_num = format(fecha, "%W"),

         dia_semana = wday(fecha)) %>% 

  filter(dia_semana==5, anio==2016|2017, mes=="enero") %>% 

  group_by(anio,dia_semana_num) %>% 

  summarise(suma_venta = sum(ventas),

            fecha = max(fecha))

# A tibble: 22 x 4

# Groups:   anio [5]

   anio  dia_semana_num suma_venta fecha              

   <chr> <chr>               <dbl> <dttm>             

 1 2016  01                   2.88 2016-01-07 00:00:00

 2 2016  02                   7.51 2016-01-14 00:00:00

 3 2016  03                   9.52 2016-01-21 00:00:00

 4 2016  04                  11.4  2016-01-28 00:00:00

 5 2017  01                 150.   2017-01-05 00:00:00

 6 2017  02                 153.   2017-01-12 00:00:00

 7 2017  03                 155.   2017-01-19 00:00:00

 8 2017  04                 158.   2017-01-26 00:00:00

 9 2018  01                 304.   2018-01-04 00:00:00

10 2018  02                 307.   2018-01-11 00:00:00

# ... with 12 more rows

Esta identificación de datos, nos permite también obtener estadísticos y operaciones semanales en función de días de semana específicos. Por ejemplo, podemos obtener la mediana de ventas de los viernes de junio y julio 2016, utilizando la función filter:

data_simulada_venta  %>%

  mutate(anio = format(fecha, "%Y"),

         mes = format(fecha, "%B"),

         dia_semana_num = format(fecha, "%W"),

         dia_semana = wday(fecha)) %>% 

  filter(dia_semana==5, anio==2016, mes %in% c("junio", "julio")) %>% 

  group_by(anio,dia_semana_num) %>% 

  summarise(suma_venta = median(ventas),

            fecha = max(fecha))

`summarise()` regrouping output by 'anio' (override with `.groups` argument)

# A tibble: 9 x 4

# Groups:   anio [1]

  anio  dia_semana_num suma_venta fecha              

  <chr> <chr>               <dbl> <dttm>             

1 2016  22                   61.0 2016-06-02 00:00:00

2 2016  23                   64.3 2016-06-09 00:00:00

3 2016  24                   66.4 2016-06-16 00:00:00

4 2016  25                   68.7 2016-06-23 00:00:00

5 2016  26                   72.9 2016-06-30 00:00:00

6 2016  27                   75.5 2016-07-07 00:00:00

7 2016  28                   78.9 2016-07-14 00:00:00

8 2016  29                   82.7 2016-07-21 00:00:00

9 2016  30                   85.9 2016-07-28 00:00:00

Una vez obtenidos los datos semanales, podemos utilizar los paquetes de pronóstico de R (forecast y tseries) para realizar proyecciones sea de las series semanales o diarias:

a.      Proyección semanal.

ts(data_semanal$suma_venta) %>%

  auto.arima() %>% 

  forecast(h=5)

Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95

230       3369.852 3031.907 3707.797 2853.010 3886.694

231       3755.108 3375.773 4134.443 3174.965 4335.251

232       3727.434 3339.778 4115.089 3134.565 4320.302

233       3685.478 3287.301 4083.655 3076.519 4294.437

234       3653.735 3241.844 4065.626 3023.802 4283.668

b.     Proyección diaria (usando los datos de la semana completa).

ts(data_simulada_venta$ventas) %>%

  auto.arima() %>% 

  forecast(h=5)

     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95

1357       674.7300 674.3008 675.1591 674.0736 675.3863

1358       675.1442 674.4392 675.8492 674.0659 676.2224

1359       675.5323 674.4952 676.5694 673.9462 677.1184

1360       675.9322 674.4860 677.3784 673.7204 678.1440

1361       676.3322 674.4527 678.2117 673.4578 679.20

c.      Proyección diaria (usando solo los lunes para proyectar las ventas de los lunes únicamente).

data_lunes <- data_simulada_venta %>%

  mutate(dia_semana_nombre = wday(fecha, label=T),

         dia_semana_num = wday(fecha),

         anio = format(fecha, "%Y"),

         mes = format(fecha, "%B")) %>% 

  filter(dia_semana_num==2) 

auto.arima(data_lunes$ventas) %>% 

  forecast(h=5)

  Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95

227       676.2422 675.3601 677.1243 674.8932 677.5913

228       679.2780 678.0324 680.5236 677.3730 681.1829

229       682.2478 680.7087 683.7869 679.8940 684.6016

230       685.3183 683.4556 687.1810 682.4695 688.1671

231       688.3682 686.2326 690.5037 685.1021 691.6342