02-principios-ejemplos.Rmd

# Principios y análisis exploratorio



## Algunos conceptos básicos {-}

Por ejemplo, los siguientes datos fueron registrados en un restaurante durante cuatro días consecutivos:

### Ejemplo {-}

```{r, echo = FALSE}
library(tidyverse)
library(ggrepel)
library(knitr)
source("./R/funciones_auxiliares.R")
propinas <- reshape2::tips %>% 
  rename(cuenta_total = total_bill, propina = tip, sexo = sex, fumador = smoker,
         dia = day, momento = time, num_personas = size) %>% 
  mutate(sexo = recode(sexo, Female = "Mujer", Male = "Hombre"), 
         fumador = recode(fumador, No = "No", Si = "Si"),
         dia = recode(dia, Sun = "Dom", Sat = "Sab", Thur = "Jue", Fri = "Vie"),
         momento = recode(momento, Dinner = "Cena", Lunch = "Comida")) %>% 
  select(-sexo) %>% 
  mutate(dia  = fct_relevel(dia, c("Jue", "Vie", "Sab", "Dom")))
```

```{r}
sample_n(propinas, 10) %>% formatear_tabla()
```

Aquí la unidad de observación es una cuenta particular. Tenemos tres mediciones numéricas
de cada cuenta: cúanto fue la cuenta total, la propina, y el número de personas asociadas a la cuenta.
Los datos están separados según se fumó o no en la mesa, y temporalmente en dos partes: el día
(Jueves, Viernes, Sábado o Domingo), cada uno separado por Cena y Comida. 

El primer tipo de comparaciones que nos interesa hacer es para una medición: ¿Varían mucho o poco los datos 
de un tipo  de medición? ¿Cuáles son valores típicos o centrales? ¿Existen valores atípicos?

Supongamos entonces que consideramos simplemente la variable de *cuenta_total*. Podemos comenzar
por **ordenar los datos**, y ver cuáles datos están en los extremos y cuáles están en los lugares centrales:

```{r}
propinas <- propinas %>% 
  mutate(orden_cuenta = rank(cuenta_total, ties.method = "first"), 
         f = (orden_cuenta - 0.5) / n()) 
cuenta <- propinas %>% select(orden_cuenta, f, cuenta_total) %>% arrange(f)
bind_rows(head(cuenta), tail(cuenta)) %>% formatear_tabla()
```


y graficamos los datos en orden, interpolando valores consecutivos.

```{r, fig.width = 7, fig.height = 4, echo = FALSE}
g_orden <- ggplot(cuenta, aes(x = orden_cuenta, y = cuenta_total)) + geom_point(colour = "red", alpha = 0.5) + 
  labs(subtitle = "Cuenta total") 
g_cuantiles <- ggplot(cuenta, aes(x = f, y = cuenta_total)) + geom_point(colour = "red", alpha = 0.5) + geom_line() +
  labs(subtitle = "") +
  scale_x_continuous(breaks = seq(0, 1, 0.1))
gridExtra::grid.arrange(g_orden, g_cuantiles, ncol = 2)
```

A esta función le llamamos la **función de cuantiles** para la variable cuenta total. Nos
sirve para comparar directamente los distintos valores que observamos los datos
según el orden que ocupan. 

**Dispersión y Valores centrales**

- El *rango* de datos va de unos 3 dólares hasta 50 dólares
- Los **valores centrales** (del cuantil 0.25 al 0.75, por ejemplo), están
entre unos 13 y 25 dólares
- Podemos usar el cuantil 0.5 (**mediana**) para dar un valor *central* de esta distribución,
que está alrededor de 18 dólares.

Y podemos dar resúmenes más refinados si es necesario

- El cuantil 0.95 es de unos 35 dólares - sólo 5\% de las cuentas son de más de 35 dólares
- El cuantil 0.05 es de unos 8 dólares - sólo 5\% de las cuentas son de 8 dólares o menos.

Finalmente, la forma de la gráfica se interpreta usando su pendientes, haciendo comparaciones
de diferentes partes de la gráfica:

- La distribución de valores tiene asimetría: el 10\% de las cuentas más altas 
tiene considerablemente más dispersión que el 10\% de las cuentas más bajas. 

- Entre los cuantiles 0.2 y 0.5 es donde existe *mayor* densidad de datos: la pendiente
es baja, lo que significa que al avanzar en los cuantiles, los valores observados no cambian mucho.

- Cuando la pendiente es alta, quiere decir que los datos tienen más dispersión local o están más separados.

En algunos casos, es más natural hacer un *histograma*, donde dividimos el rango de la variable
en cubetas o intervalos (en este caso de igual longitud), y graficamos cuántos datos caen en cada
cubeta:

```{r, fig.width = 10, fig.height = 4, echo = FALSE}
variable <- quo(cuenta_total)
binwidth_min = 1
g_1 <- ggplot(propinas, aes(x = !!variable)) + geom_histogram(binwidth = binwidth_min) 
g_2 <- ggplot(propinas, aes(x = !!variable)) + geom_histogram(binwidth = binwidth_min * 2) 
g_3 <- ggplot(propinas, aes(x = !!variable)) + geom_histogram(binwidth = binwidth_min * 5) 
gridExtra::grid.arrange(g_1, g_2, g_3, ncol = 3)
```

Es una gráfica más popular, pero perdemos cierto nivel de detalle, y distintas particiones 
resaltan distintos aspectos de los datos.

Finalmente, una gráfica más compacta que resume la gráfica de cuantiles o el histograma
es el diagrama de caja y brazos. Mostramos dos versiones, la clásica de Tukey (T) y otra versión
menos común de Spear/Tufte (ST):

```{r, fig.width = 8, fig.height = 4, echo = FALSE}
library(ggthemes)
cuartiles <- quantile(cuenta$cuenta_total)
g_1 <- ggplot(cuenta, aes(x = f, y = cuenta_total)) + 
  labs(subtitle = "Gráfica de cuantiles: Cuenta total") +
  geom_hline(yintercept = cuartiles[2], colour = "gray") + 
  geom_hline(yintercept = cuartiles[3], colour = "gray") +
  geom_hline(yintercept = cuartiles[4], colour = "gray") +
  geom_point(alpha = 0.5) + geom_line() 
g_2 <- ggplot(cuenta, aes(x = factor("ST", levels =c("ST")), y = cuenta_total)) + 
  geom_tufteboxplot() +
  #geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
  labs(subtitle = " ") +  xlab("") + ylab("")
g_3 <- ggplot(cuenta, aes(x = factor("T"), y = cuenta_total)) + geom_boxplot() +
  labs(subtitle = " ") +  xlab("") + ylab("")
g_4 <- ggplot(cuenta, aes(x = factor("P"), y = cuenta_total)) + geom_jitter(height = 0, width =0.2, alpha = 0.5) +
  labs(subtitle = " ") +  xlab("") + ylab("")
gridExtra::grid.arrange(g_1, g_2, g_3, g_4, nrow = 1, widths = c(8, 2, 2, 2)) 
```


**Ventajas en el análisis inicial**

En un principio del análisis, estos resúmenes
(cuantiles) pueden ser más útiles que utilizar medias y varianzas, por ejemplo. La razón es
que los cuantiles:

- Son cantidades más fácilmente interpretables
- Los cuantiles centrales son más resistentes a valores atípicos que medias o varianzas
- Sin embargo, permite identificar valores extremos
- Es fácil comparar cuantiles de distintos bonches de datos



### Ejercicio {-}

¿Cómo se ve la gráfica de cuantiles de las propinas? ¿Cómo crees que esta gráfica se
compara con distintos histogramas?

```{r, fig.width = 4, fig.height = 3}
g_1 <- ggplot(propinas, aes(sample = propina)) + 
  geom_qq(distribution = stats::qunif) + xlab("f") + ylab("propina")
g_1
```

**Observación**. Cuando hay datos repetidos, los cuantiles tienen que interpretarse como sigue: 
el cuantil-$f$ con valor $q$ satisface que existe una proporción aproximada $f$ de los datos que están en el valor
$q$ o por debajo de éste, pero no necesariamente exactamente una proporción $f$ de los datos estan en $q$ o por debajo.

## Ejemplo: Precio de venta de casas {-}

Consideramos [datos de precios de ventas de la ciudad de Ames, Iowa](https://www.kaggle.com/prevek18/ames-housing-dataset).
Nos interesa entender la variación del precio de las casas.

```{r, echo = FALSE, message = FALSE, warning = FALSE}
source("R/casas_preprocesamiento.R")
set.seed(21)
casas_completo <- casas
casas <- casas_completo %>% sample_frac(0.9)
casas_holdout <- casas_completo %>% anti_join(casas)
nombres <- casas %>% group_by(nombre_zona) %>% tally %>% pull(nombre_zona)
zonas <- nombres[table(casas$nombre_zona) > 30]
casas <- casas %>% mutate(nombre_zona = fct_reorder(nombre_zona, precio_miles)) %>% 
  filter(nombre_zona %in% zonas) %>% 
  mutate(precio_m2 = precio_m2_miles * 1000)
```

Calculamos primeros unos cuantiles de los precios de las casas:

```{r}
quantile(casas %>% pull(precio_miles)) 
```

Una primera comparación que podemos hacer es considerar las distintas zonas de la ciudad.
Podemos usar diagramas de caja y brazos para **comparar** precios en distintas zonas
de la ciudad:

```{r}
ggplot(casas, aes(x = nombre_zona, y = precio_miles)) + geom_boxplot() + coord_flip()
```

La primera pregunta que nos hacemos de esta comparación es cómo pueden variar características
de las casas dentro de cada zona. En primer lugar, podemos considerar el área de las casas. En lugar
de graficar el precio, graficamos el precio por metro cuadrado, por ejemplo:

```{r, echo = FALSE}
casas <- casas %>% mutate(nombre_zona = fct_reorder(nombre_zona, precio_m2_miles)) %>% 
  filter(nombre_zona %in% zonas)
```

```{r}
ggplot(casas, aes(x = nombre_zona, y = precio_m2)) + geom_boxplot() + coord_flip()
```


Podemos cuantificar la variación que observamos de zona a zona y la variación que hay dentro de zonas. La 
variación que vemos entre las medianas de la zona es:

```{r}
casas %>% group_by(nombre_zona) %>% summarise(mediana_zona = median(precio_m2)) %>% 
  pull(mediana_zona) %>% quantile %>% round
```

Y las variaciones con respecto a las medianas **dentro** de cada zona, agrupadas, se resume como:

```{r}
quantile(casas %>% group_by(nombre_zona) %>% 
           mutate(residual = precio_m2 - median(precio_m2)) %>% pull(residual)) %>% round
```

Nótese que este último paso tiene sentido pues la variación dentro de las zonas, en términos de precio por metro
cuadrado, es similar. Esto no lo podríamos hacer de manera efectiva si hubiéramos usado el precio de las casas sin
ajustar por su tamaño.

Podemos resumir este primer análisis con un par de gráficas de cuantiles (@ClevelandVis):

```{r, fig.width = 5, fig.height = 3}
mediana <- median(casas$precio_m2)
resumen <- casas %>% 
  group_by(nombre_zona) %>% mutate(mediana_zona = median(precio_m2)) %>% 
  mutate(residual = precio_m2 - mediana_zona) %>% ungroup %>% 
  mutate(mediana_zona = mediana_zona - mediana) %>% 
  select(nombre_zona, mediana_zona, residual) %>% 
  gather(tipo, valor, mediana_zona:residual) 
ggplot(resumen, aes(sample = valor)) + geom_qq(distribution = stats::qunif) + 
  facet_wrap(~ tipo) + ylab("Precio por m2") + xlab("f") +
  labs(subtitle = "Precio por m2 por zona", caption = paste0("Mediana total de ", round(mediana)))
```

Y vemos que la mayor parte de la variación del precio por metro cuadrado ocurre dentro de cada
zona, una vez que controlamos por el tamaño de las casas. La variación dentro de cada zona
es aproximadamente simétrica, aunque la cola derecha es ligeramente más larga con algunos valores
extremos.

Podemos seguir con otro indicador importante: la calificación de calidad de los terminados
de las casas. Como primer intento podríamos hacer:

```{r, echo = FALSE}
ggplot(casas %>% mutate(ind_calidad = cut(calidad_gral, c(0, 5, 8,10))), aes(y = precio_m2, x = nombre_zona, 
        colour=factor(ind_calidad))) + 
  geom_hline(yintercept = c(1000, 2000), colour = "gray") +
  geom_jitter(width = 0.2, height = 0, alpha = 0.5) +
  coord_flip() + scale_colour_colorblind() +
  facet_wrap(~ind_calidad)

```

Lo que indica que las calificaciones de calidad están distribuidas de manera
muy distinta a lo largo de las zonas, y que probablemente no va ser simple
desentrañar qué variación del precio se debe a la zona y cuál se debe a la calidad.



## Enlace {-}

Consideremos la prueba Enlace (2011) de matemáticas para primarias. Una primera pregunta 
que alguien podría hacerse es:  ¿qué escuelas son mejores, las privadas o las públicas? 

```{r, message = FALSE, echo = FALSE, include=FALSE, eval = FALSE}
#col_spec <- cols_only(X2 = col_character(), X3 = col_character(),
#                      X6 = col_character(), X24 = col_double(),
#                      X82 = col_integer(),
#                      X84 = col_character())
#enlace_1 <- read_csv("./datos/enlace/escuelas_enlace_nacional_primaria_1.csv", skip=3, 
#                     col_names = FALSE, guess_max = 100000, col_type = col_spec,
#                     locale = locale(encoding = "ISO-8859-1"))
#enlace_2 <- read_csv("./datos/enlace/escuelas_enlace_nacional_primaria_2.csv", skip=3, 
#                     col_names = FALSE, guess_max = 100000, col_type = col_spec,
#                     locale = locale(encoding = "ISO-8859-1"))
#enlace <- bind_rows(enlace_1, enlace_2)
#names(enlace) <- c("estado", "clave", "tipo", "mate_6",  "num_evaluados_total", "marginacion")
#table(enlace$tipo)
#enlace <- enlace %>% mutate(tipo = ifelse(tipo == "INDêGENA", "INDÍGENA", tipo))
#write_csv(enlace, "datos/enlace.csv")
enlace <- read_csv("datos/enlace.csv")
```

Con estos datos a la mano, podemos hacer unos primeros resúmenes de los datos. El rango
de la calificación de matemáticas para un alumno es de 0-800, y aproximadamente
la mitad de los alumnos califica en el rengo de 450 a 550. Vemos dispersión
considerable en las calificaciones de las escuelas, y diferencias considerables 
entre tipo de escuelas:

```{r, message = FALSE, echo = FALSE}
enlace <- read_csv("datos/enlace.csv")
enlace <- enlace %>%  filter(num_evaluados_total > 0, mate_6 > 0) %>%
    mutate(marginacion = fct_reorder(marginacion, mate_6, median)) %>% 
    mutate(tipo = recode(tipo, `INDÍGENA`="Indígena/Conafe", CONAFE="Indígena/Conafe", GENERAL="General", PARTICULAR="Particular")) %>% 
    mutate(tipo = fct_reorder(tipo, mate_6, mean)) 
```

```{r, message = FALSE}
enlace_tbl <- enlace %>% group_by(tipo) %>% 
    summarise(n_escuelas = n(),
              cuantiles = list(cuantil(mate_6, c(0.05, 0.25, 0.5, 0.75, 0.95)))) %>% 
    unnest(cols = cuantiles) %>% mutate(valor = round(valor)) 
enlace_tbl %>% spread(cuantil, valor) %>% formatear_tabla()
```

Podemos graficar de varias maneras, por ejemplo (en la gráfica 1 usamos cuantiles 0.05, 0.25, 0.5, 0.75 y 0.95) :


```{r, fig.width = 10, fig.height = 6, echo = FALSE, message = FALSE}
g_medianas <- ggplot(enlace_tbl %>% filter(cuantil == 0.50), aes(x = tipo, y = valor)) +
    geom_point(colour = "red") + ylim(c(150,880)) + labs(subtitle = "Gráfica 1")
g_80 <- ggplot(enlace_tbl  %>% spread(cuantil,valor), 
                                     aes(x = tipo, y = `0.5`)) +
    geom_linerange(aes(ymin= `0.05`, ymax = `0.95`), colour = "gray40") +
    geom_point(colour = "red", size = 3) + ylim(c(150,880)) + 
    labs(subtitle = "Gráfica 1") +
    ylab("Promedios Matemáticas")
    
g_80_p <- ggplot(enlace_tbl  %>% spread(cuantil,valor), aes(x = tipo, y = `0.5`)) +
    geom_linerange(aes(ymin= `0.05`, ymax = `0.95`), colour = "gray40") +
    geom_linerange(aes(ymin= `0.25`, ymax = `0.75`), size = 2, colour = "white") +
    geom_point(colour = "red", size = 3) +
     ylim(c(150,880)) + labs(subtitle = "Gráfica 1")+
    ylab("Promedios Matemáticas")

g_boxplot <- ggplot(enlace  , aes(x = tipo, y = mate_6)) + 
    geom_boxplot(outlier.size = 0.7) +
    labs(subtitle = "Gráfica 2")+ ylim(c(150,930))+
    ylab("Promedios Matemáticas")

g_cuantil <- ggplot(enlace, aes(sample = mate_6, colour = tipo)) +
  geom_qq(distribution = stats::qunif, size = 1) + ylab("Promedios Matemáticas") +
  xlab("orden") + labs(subtitle = "Gráfica 3")

gridExtra::grid.arrange(g_80_p, g_boxplot, g_cuantil, nrow = 2)
```

En términos de comparaciones, podemos discutir qué tan apropiada es cada gráfica. Graficar más
cuantiles es más útil para hacer comparaciones.  Por ejemplo, en la Gráfica 1 podemos ver que la
mediana de las escuelas generales está cercano al cuantil 5\% de las escuelas particulares. Por otro
lado, el diagrama de caja y brazos muestra también valores "atípicos". 

La diferencia es considerable entre tipos de escuela, pero antes de contestar prematuramente
la pregunta: ¿cuáles son las mejores escuelas?
busquemos mejorar la interpretabilidad de nuestras comparaciones usando los principios 2 y 3. Podemos comenzar
por agregar, por ejemplo, el nivel del marginación del municipio donde se encuentra la escuela.

```{r, echo = FALSE}
enlace_tbl_marg <- enlace %>% 
    group_by(tipo, marginacion) %>% 
    summarise(n_alumnos = sum(num_evaluados_total),
              cuantiles = list(cuantil(mate_6, c(0.05, 0.25, 0.5, 0.75, 0.95)))) %>% 
    unnest(cols = cuantiles) %>% mutate(valor = round(valor)) %>% 
    filter( n_alumnos > 20)
```

Podemos graficar de esta manera, usándo páneles (pequeños múltiplos útiles para hacer comparaciones):

```{r, fig.width = 10, fig.height = 4, echo = FALSE}
g_80_p <- ggplot(enlace_tbl_marg  %>% spread(cuantil, valor), 
                                     aes(x = marginacion, y = `0.5`)) +
    geom_linerange(aes(ymin= `0.05`, ymax = `0.95`), colour = "gray40") +
    geom_linerange(aes(ymin= `0.25`, ymax = `0.75`), size = 2, colour = "white") +
    geom_point(colour = "red", aes(size = log10(n_alumnos/1000))) +
    ylab("Promedios Matemáticas") + facet_wrap(~tipo, nrow = 1) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    scale_size_continuous(name = "Miles de \nAlumnos",
                          breaks = c(-0.3, 0, 1, 2, 2.7),
                          labels = c(0.5, 1, 10, 100, 500)) 
g_80_p
```



Esta gráfica pone en contexto la pregunta inicial, y señala la dificultad 
de contestarla. En primer lugar,

- Señala que la pregunta no sólo debe concentarse en el tipo de "sistema": pública, privada, etc. Las
escuelas públicas en zonas de marginación baja no tienen una distribución de calificaciones muy distinta
a las privadas en zonas de marginación alta
- Quiere decir que el contexto de la escuela es importante.
- Y si eso es importante, podríamos entonces pensar, por ejemplo, que factores como el entorno familiar
de los estudiantes (por ejemplo: cosas tan extremas como si tienen que trabajar o no) puede resultar en comparaciones
que favorecen a las escuelas privadas - cuando quizá en parte no se debe a que tengan un sistema especialmente
bueno.
- Si esto es cierto, entonces la pregunta inicial es demasiado vaga y mal planteada: quizá deberíamos intentar
entender cuánto "aporta" cada escuela a cada estudiante, como medida de qué tan buena es cada escuela.



## Estados y calificaciones en SAT {-}

¿Cómo se relaciona el gasto por alumno, a nivel estatal, 
con sus resultados académicos? Hay trabajo
considerable en definir estos términos, pero supongamos que tenemos el
[siguiente conjunto de datos](http://jse.amstat.org/datasets/sat.txt) (@Guber2011, @sleuth), que son
datos oficiales agregados por estado de Estados Unidos. Tenemos las variables
*sat*, por ejemplo, que es la calificación promedio de los alumnos en cada estado
(para 1997), y la variable *expend*, que es el gasto en miles de dólares
por estudiante en (1994-1995), además de algunas otras variables. 


```{r, message = FALSE}
sat <- read_csv("datos/sat.csv")
sat_tbl <- sat %>% select(state, expend, sat) %>% 
    gather(variable, valor, expend:sat) %>% 
    group_by(variable) %>% 
    summarise(cuantiles = list(cuantil(valor))) %>% 
    unnest(cols = c(cuantiles)) %>% 
    mutate(valor = round(valor, 1)) %>% 
    spread(cuantil, valor)
sat_tbl %>% formatear_tabla
```


Esta variación considerable es considerable para promedios del SAT: 
el percentil 75 es alrededor de 1050 puntos, mientras que el percentil 33 corresponde a alrededor de 800, 
e igualmente hay diferencias considerables de gasto por alumno (miles de dólares) a lo largo 
de los estados.

Ahora hacemos nuestro primer ejercico de comparación: ¿Cómo se ven las
calificaciones para estados en distintos niveles de gasto? Podemos
usar una gráfica de dispersión:


```{r}
 ggplot(sat, aes(x = expend, y = sat, label = state)) + 
  geom_point(colour = "red", size = 2) + geom_text_repel(colour = "gray50") +
  xlab("Gasto por alumno (miles de dólares)") +
  ylab("Calificación promedio en SAT")
```

Estas comparaciones no son de muy alta calidad, solo estamos usando 2 variables (pocas),
y no hay mucho que podamos decir en cuanto explicación.

**Las unidades que estamos comparando pueden diferir fuertemente en otras 
dimensiones importantes, lo cual hace interpretar la gráfica muy difícil**

Sabemos que es posible que el IQ difiera en los estados, pero no como producir
diferencias de este tipo. Sin embargo, descubrimos que existe una variable adicional, 
que es el porcentaje de alumnos de cada estado
que toma el SAT. Podemos agregar como sigue:

```{r}
 ggplot(sat, aes(x = expend, y = math, label=state, colour = frac)) + 
  geom_point() + geom_text_repel() +
  xlab("Gasto por alumno (miles de dólares)") +
  ylab("Calificación en matemáticas") 
```


Y vemos entonces por qué nuestra comparación inicial es relativamente pobre:
los estados con mejores resultados promedio en el SAT son aquellos donde una
fracción relativamente baja de los estudiantes toma el examen. La diferencia
es considerable.

En este punto podemos hacer varias cosas. Una primera idea es intentar comparar
estados más similares en cuanto a la población de alumnos que asiste. Podríamos hacer
grupos como sigue:

```{r, echo = FALSE, fig.width = 6, fig.height = 3}
set.seed(991)
sat$clase <- kmeans(sat %>% select(frac), 
                    centers = 4,  nstart = 100, iter.max = 100)$cluster
sat <- sat %>% group_by(clase) %>% mutate(clase_media = round(mean(frac))) %>% 
  ungroup %>% mutate(clase_media = factor(clase_media))
sat <- sat %>% mutate(rank_p = rank(frac, ties= "first") / length(frac))
ggplot(sat, aes(x = rank_p, y = frac, label = state, 
                colour = clase_media)) +
  geom_point(size = 2) + 
  #geom_text_repel(colour = "gray80", size = 3) +
  paleta
```

Estos resultados indican que es más probable que buenos alumnos decidan hacer
el SAT - y esto ocurre de manera diferente en cada estado (en algunos estados
era más común otro examen, el ACT). 

Si hacemos clusters de estados
según el % de alumnos, empezamos a ver otra historia. Ajustamos rectas
de mínimos cuadrados como referencia:

```{r, echo = FALSE}
ggplot(sat, aes(x = expend, y = math, label=state, colour = clase_media)) + 
  geom_point(size = 3) + 
  geom_smooth(method = "lm", se = FALSE) +
  xlab("Gasto por alumno (miles)") +
  ylab("Calificación en matemáticas") + paleta +
  geom_text_repel(colour = "gray70") 
```



Sin embargo, el resultado puede variar considerablemente si categorizamos de distintas maneras.

Otra idea es más útil es *factorizar* o *controlar*
el efecto de el % de alumnos que toma el examen, y hacer una comparación de la variación
restante. Para hacer esto, podemos examinar primero la relación entre la fracción de
alumnos que toma el examen del SAT y las calificaciones:

```{r, fig.width = 5, fig.height = 3, echo = FALSE}
ggplot(sat, aes(x = frac, y = sat, label = state)) + geom_point() +
    geom_smooth(method = "loess", span = 0.5, se = FALSE, method.args = list(degree = 1)) +
    annotate("text", x = 40, y = 1050, label = "Desempeño por \n encima del esperado", colour = "gray30") +
    annotate("text", x = 15, y = 900, label = "Desempeño por \n abajo del esperado", colour  = "gray30") 

```

Una buena parte de variación de la calificación del sat puede explicarse con las la proporción 
de estudiantes que toma el examen. Ahora podemos graficar los residuales
obtenidos de este ajuste contra el gasto por alumno (a los residuales le
sumamos la mediana de las calificaciones del SAT). A estos residuales
le podemos llamar *calificación SAT ajustada por porcentaje de alumnos que
toma el examen):



```{r, echo = FALSE, message = FALSE}	
suavizador <- loess(sat ~ frac, sat, degree = 1, span = 0.5)	
sat <- sat %>% mutate(residual = residuals(suavizador), ajustado = fitted(suavizador))	
mediana <- median(fitted(suavizador))	
# ggplot(sat %>% select(ajustado, residual, state) %>% gather(tipo, valor, residual:ajustado),	
#   aes(sample = valor - (tipo=="ajustado")*mediana)) + geom_qq(distribution = stats::qunif) + 	
#   facet_wrap(~tipo)	
```


```{r}
 ggplot(sat, aes(x = expend, y = residual + mediana, label = state)) + 
  geom_point(colour = "red", size = 2) + geom_text_repel(colour = "gray50") +
  xlab("Gasto por alumno (miles de dólares)") +
  ylab("Ajustado SAT (vs fracción de alumnos)")
```

La indicación que obtenemos es diferente: ahora parece haber cierta 
evidencia de  que estados con más gasto
tienden a tener mejores resultados, aunque  a total, la variación
del SAT ajustado es considerablemente menor que la del SAT no ajustado.

Nótese que esto nos muestra que nuestra pregunta original tiene una respuesta
simple pero poco útil. con este análisis podemos pensar en mejores preguntas, pero todas
son mucho más difíciles de contestar!

### Ejercicio: suavizadores {-}

Puedes practicar y entender cómo funcionan los suavidaores con **esta app**, cuya idea
es tomada de [este libro](https://rafalab.github.io/dsbook/smoothing.html#local-weighted-regression-loess). El suavizamiento
**loess** (@ClevelandVis) es diferente de promedios móviles centrados usuales, y a veces ayuda a obtener mejores
suavizados (cuando hay estructura lineal en el problema):

```{r, fig.width = 10, fig.height = 6, echo = FALSE}

w_loess <- function(x){
  ifelse(abs(x) < 1, (1 - abs(x)^3)^3, 0.0)
}
lm_ponderado <- function(...) {
  geom_smooth(method = "lm", se = FALSE, ...)
}

crear_local <- function(frac_val = 10, alpha = 0.5){
  num_datos = floor(alpha * nrow(sat)) + 1
  sat_g <- sat %>% 
    mutate(activos = rank(abs(frac - frac_val)) < num_datos,
         max_dist = max(activos*abs(frac - frac_val)),
         peso = w_loess((frac - frac_val)/max_dist))
ggplot(sat_g, aes(x = frac, y = sat, label = state, size = peso)) + 
  geom_point(colour = "gray") +
  scale_size(range = c(0.5, 4)) +
    geom_smooth(method = "loess", method.args = list(degree = 1), span = alpha, se = FALSE,
              size = 1, colour = "red") +
  lm_ponderado(data = sat_g %>% filter(peso > 0), aes(weight = peso^2)) +
  geom_vline(xintercept = frac_val, colour = "gray") +
  theme(legend.position = "none") 

}
grafs <- map(c(10, 20, 30, 40, 50, 60), ~ crear_local(.x, alpha = 0.6))
gridExtra::grid.arrange(grobs = grafs,ncol = 3)
```


Promedios móviles (ponderados por distancia) dan resultados más ruidosos si queremos
capturar la estructura lineal local:

```{r, fig.width = 10, fig.height = 6, echo = FALSE, echo = FALSE}

w_loess <- function(x){
  ifelse(abs(x) < 1, (1 - abs(x)^3)^3, 0.0)
}
lm_ponderado <- function(...) {
  geom_smooth(method = "lm", formula = y ~ 1, se = FALSE, ...)
}

crear_local <- function(frac_val, alpha = 0.5){
  num_datos = floor(alpha * nrow(sat)) 
  sat_g <- sat %>% 
    mutate(activos = rank(abs(frac - frac_val)) < num_datos,
         max_dist = max(activos*abs(frac - frac_val)),
         peso = w_loess((frac - frac_val)/max_dist))
ggplot(sat_g, aes(x = frac, y = sat, label = state, size = peso)) + 
  geom_point(colour = "gray") +
  scale_size(range = c(0.5, 4)) +
    geom_smooth(method = "loess", method.args = list(degree = 0), 
                span = alpha, se = FALSE,
              size = 1, colour = "red") +
  lm_ponderado(data = sat_g %>% filter(peso > 0), aes(weight = peso)) +
  geom_vline(xintercept = frac_val, colour = "gray") +
  theme(legend.position = "none") 

}
grafs <- map(c(10, 20, 30, 40, 50, 60), ~ crear_local(.x, alpha = 0.3))
gridExtra::grid.arrange(grobs = grafs,ncol = 3)
```



## Ajuste y residuales {-}

Esta técnica de análisis en la que controlamos por variables para "escarbar en los datos",
se puede resumir como sigue. Intentamos ajustar patrones conocidos, "quitarlos de
los datos", y examinar la variación restante (@ClevelandVis):

$$ datos = ajuste + residual $$

En un principio, ajuste y residuales requieren de nuestra atención. El objetivo de este proceso
es hacer más datos, no menos, con la idea de que hacer comparaciones más finas que las que los
datos sin procesar nos dan. Es una técnica que repetiremos varias veces.



## Comparaciones multiplicativas {-}

Hay dos maneras de hacer comparaciones:

- En escala **aditivas**: por ejemplo, el grupo A tiene 50 personas más que el grupo B
- En escala **multiplicativas**: el grupo A tiene 25\% más personas que el grupo B

¿Cuáles son más apropiadas?

Distintos aspectos del análisis aparecen dependiendo de qué tipo de comparaciones queramos
hacer, pero muy frecuentemente una tipo comparación es preferible al otro. Por ejemplo,
si queremos comparar incrementos de precios a lo largo de distintos tipos de productos, es
natural usar escalas multiplicativas. 

La elección de la comparación apropiada mejora la calidad de las comparaciones, y
simplifica el análisis y su interpretación. 

En el siguiente ejemplo, consideramos cómo varía el precio de un conjunto de casas
dependiendo de la calificación de calidad:



```{r, echo = FALSE, message =FALSE}
casas_ejemplo <- casas %>% filter(calidad_gral > 3) %>% 
           mutate(nombre_zona = fct_reorder(nombre_zona, calidad_gral))
ggplot(casas_ejemplo, aes(x = factor(calidad_gral), y = precio_miles)) + 
    geom_boxplot(outlier.alpha = 0) +
    geom_jitter(width =0.1, height = 0, alpha = 0.2) 
```

En esta gráfica podríamos hacer comparaciones entre precios de casas de manera
aditiva:

- Por ejemplo, cuando aumentamos de calificación 7 a 8, la mediana de precio
sube en unos 75 dólares por metro cuadrado. Un aumento de 4 a 5 sin embargo,
representa alrededor de  20 dólares adicionales.
- A partir de la mediana, la variación intercuartil de los precios cuando
la calidad es 4 es de más menos  10 dólares. En calidad 10, la variación es
de alrededor de más menos 40 dólares
- Podemos seguir así sucesivamente, pero nuestra descripción **aditiva** es complicada,
y no muestra claramente alguna estructura consistente.

Sin embargo, una comparación **multiplicativa** puede ser más adecuada y simplifica la descripción
de los datos. Por ejemplo, una variación de 10 mil dólares en casas de calidad 10 no es 
muy considerable (con rango de 200 a 600 mil dólares), pero es muy importante para una casa de calidad 4.

Podemos transformar nuestros datos para trabajar naturalmente con ellos en una escala que privilegie comparaciones multiplicativas. Una función apropiada es el logaritmo
natural, pues es una función que cumple:

$$\log(y) - \log(x) = \log(y/x) = \log \left(1 + \frac{y - x}{x}\right) \approx \frac{y - x}{x}$$

y esta aproximación es buena cuando $\left |\frac{y-x}{x} \right | < 0.25$ (y muy buena si es menor a 0.1). Es decir, el logartimo convierte comparaciones relativas
en comparaciones aditivas. Es más simple trabajar con comparaciones aditivas:

Veamos entonces cómo se ven los datos de casas en escala logarítmica. A la izquierda
mostramos una gráfica más apropiada para análisis, y a la derecha una que quizá prefereríamos
usar para presentación

```{r, fig.width = 9, fig.height = 3, echo = FALSE}
g_1 <- ggplot(casas_ejemplo, aes(x = factor(calidad_gral), y = log(precio_miles))) +
    geom_boxplot(outlier.alpha = 0) + 
    geom_jitter(width =0.1, height = 0, alpha = 0.1)  + 
    scale_y_continuous(breaks = seq(3, 7, 0.5))
g_2 <- ggplot(casas_ejemplo, 
    aes(x = factor(calidad_gral), y = precio_miles)) +
    geom_boxplot(outlier.alpha = 0) + 
    geom_jitter(width =0.1, height = 0, alpha = 0.1) + 
    scale_y_log10(breaks = c(50, 100, 200, 400)) 
gridExtra::grid.arrange(g_1, g_2, ncol = 2)
```

Nótese que en los datos transformados:

- La variación es más similar a lo largo de cada categoría de calidad. En la escala logarítmica,
las cajas (cuartil inferior a cuartil superior) está alrededor de tamaño 0.25, y la mediana
está aproximadamente centrada en cada caja. Así que la variación entre el cuartil superior
e inferior respecto a la mediana es de 12\% aproximadamente. *Esto es aproximadamente
cierto en cada grupo de calidad*. 

- Los incrementos en precio son más o menos similares cuando la calidad aumenta en un 1 punto ( de la gráfica con unidades log, más o menos 0.25, aproximadamente 25\% de incremento en la escala original.

Esta es una comparación superior, pues descubrimos una regularidad importante en los
datos (principio 3). En este caso, lo logramos cambiando a una escala más apropiada para
el problema.

Ahora consideramos la desomposicion que hicimos:

$$ precio (log) = mediana\_calidad + residual $$

y analizamos cada parte. Podríamos calcular directamente
cambios relativos cuando aumenta un punto de calidad en las medianas:

```{r}
medianas <- casas_ejemplo %>% group_by(calidad_gral) %>% 
    summarise(mediana = median(precio_miles)) %>% 
    arrange(desc(mediana)) %>% 
    mutate(incremento_porc = (mediana / lead(mediana) - 1) %>% round(2))
medianas %>% formatear_tabla
```

Y como los residuales tienen distribución similar en cada grupo 
es razonable  agrupamos para resumir:

```{r, echo = FALSE, message = FALSE}
residuales <- casas_ejemplo %>% left_join(medianas) %>% 
    group_by(calidad_gral) %>% 
    mutate(residual = log(precio_miles) - log(mediana)) 
f <- c(0, 0.10, 0.25, 0.5, 0.75, 0.9, 1)
quantile(residuales$residual, f) %>% round(2) %>% 
  formatear_tabla()
```

Así que nuestro resumen podría ser

- Con cada punto adicional de calidad, la mediana del precio aumenta alrededor de 25\%. - Dentro de cada grupo
de calidad, el 50\% de las casas está a una diferencia máxima de +/- 12% del precio mediano de la categoría,
y 80\% de las casas alrededor de +/- 27%. 

Esta es una descripción compacta y que pone en evidencia la estructura en la relación de estas dos variables.

### Otros tipos de resumen {-}

Otros de resúmenes que se pueden hacer y son comunes son:

- La media en lugar de la mediana para resumir la ubicación general de una distribución.
- La desviación estándar en lugar del rango entre cuantiles para resumir la dispersión.

La desviación estándar se calcula como sigue: si $\bar{x}$ es la media de un 
conjunto de mediciones, su desviación estándar es 

$$\hat{\sigma} = \sqrt{\frac{1}{n}\sum_{i = 1}^n (x_i - \bar{x})^2}$$

Está son medidas más especializadas y difíciles de interpretar: por ejemplo, funcionan mejor cuando las distribuciones que queremos describir son simétricas y no tienen colas
muy largas o atípicos. Tienen muchas ventajas en ciertos casos, pero en
anáisis exporatorio tendemos a usarlas menos.


## Ejemplo: tiempos de fusión {-}

Veamos el siguiente ejemplo, que es un experimento donde se midió el tiempo
que tardan distintas personas en fusionar un [estereograma](https://en.wikipedia.org/wiki/Autostereogram#/media/File:Stereogram_Tut_Random_Dot_Shark.png) para ver una imagen 3D. (@ClevelandVis).

Existen dos condiciones: en una se dio indicaciones de qué figura tenían que
buscar (VV) y en otra no se dio esa indicación. ¿Las instrucciones verbales
ayudan a fusionar más rápido el estereograma? 

```{r, message = FALSE, echo = FALSE, fig.width =4, fig.height = 3}
fusion <- read_delim("./datos/fusion_time.txt", delim = " ", trim_ws = TRUE)
ggplot(fusion, aes(x = nv.vv, y = time)) + geom_boxplot() + geom_jitter()
```

Nótese que aquí no nos procupamos, como en los ejemplos anteriores, por diferencias
sistemáticas producidas por otros factores que los que consideramos en esta
comparación. El hecho de que se aleatorizó el tratamiento excluye esta posibiidad en gran parte.

Lo primero que observamos es que hay una variabilidad grande entre el tiempo que tardan
en fusionar: los rangos van de menos de 5 segundos hasta 20 - 40 segundos. Esto sugiere que una comparación multiplicativa puede ser más apropiada para comparar. La razón es 
que 

- No pensamos que el efecto de la información verbal, por ejemplo, cambie el tiempo de fusión en un número fijo de segudos.
- Más bien, hay personas que tienden a tardar mucho más que otras. El efecto de la instrucción verbal entonces tiene más sentido como un efecto multipicativo: quizá reduce en 20\%, o 40\%, etc el tiempo de detección. 

Decidimos entonces convertir a escala multiplicativa:

```{r, message = FALSE, echo = FALSE, fig.width =4, fig.height = 3}
fusion <- read_delim("./datos/fusion_time.txt", delim = " ", trim_ws = TRUE)
g_1 <- ggplot(fusion, aes(x = nv.vv, y = log(time))) + geom_boxplot() + geom_jitter() 
g_2 <- ggplot(fusion, aes(x = nv.vv, y = time)) + geom_boxplot() + geom_jitter() +
     scale_y_log10()
gridExtra::grid.arrange(g_1, g_2, ncol = 2)
```

De la gráfica, vemos que la distribución VV está desplazada alrededor de 0.3 unidades (en logaritmo)
hacia abajo, o en la escala original, alrededor de 35\% (exp(0.3)). 
Esto implica que los tiempos bajo VV son en promedio alrededor de 66% los tiempos bajo NV. La variación dentro de cada grupo es considerable.

Gráficas de cuantiles puedes ser muy efectivas para mostrar esta descripción:

```{r, , fig.width =4, fig.height = 3}
ggplot(fusion, aes(sample = log(time), colour = nv.vv)) + 
  geom_qq(distribution = stats::qunif)  + ylab("Log tiempo de fusión") +
  xlab("f")
```



## Cuándo la elección es menos importante {-}

En algunos casos, la diferencia entre hacer análisis con logaritmo o en unidades
originales no es muy importante. Esto sucede cuando el rango de variación está dentro de un orden de magnitud. Esto es porque cuando $x$ no varía mucho, por ejemplo,
las comparaciones aditivas $(y-x)$ y multiplicativas $(y-x)/x$ difieren en aproximadamente por un factor multiplicativo.



## Tablas de conteos {-}

Para tablas de conteos o porcentajes, las comparaciones multiplicativas tienden a ser más útiles. 
Consideremos los siguientes datos de tomadores de te (del paquete @factominer):

```{r, echo=FALSE, message=FALSE}
tea <- read_csv(("datos/tea.csv"))
te <- tea %>% select(how, price, sugar)
```

```{r}
sample_n(te, 10)
```

Nos interesa ver qué personas compran té suelto, y de qué tipo. Empezamos por ver las proporciones
que compran té según su empaque (en bolsita o suelto):

```{r}
precio <- te %>% group_by(price) %>% 
  tally() %>% mutate(prop = round(100 * n / sum(n))) %>% 
  select(-n)
tipo <- te %>% group_by(how) %>% tally() %>% 
  mutate(`%` = round(100 * n / sum(n)))
tipo %>% formatear_tabla
```

La mayor parte de las personas toma té en bolsas. Sin embargo, el tipo de té  (en términos de precio o marca) que compran es
muy distinto dependiendo de la presentación


```{r, echo = FALSE}
tipo <- tipo %>% select(how, prop_how = `%`)
tabla_2 <- te %>% group_by(how, price) %>% tally() %>% group_by(how) %>% mutate(prop = round(100 * n / sum(n))) %>% 
  select(-n) 
tabla_2 %>% spread(how, prop, fill = 0) %>% formatear_tabla()
```

Estos datos podemos examinarlos un rato y llegar a conclusiones, pero no son muy efectivo en mostrar
claramente los patrones. Tampoco es la siguiente gráfica:

```{r}
ggplot(tabla_2 %>% ungroup %>% mutate(price = fct_reorder(price, prop)),
  aes(x = price, y = prop, group = how, colour = how)) + 
  geom_point() + coord_flip() + geom_line()
```

En lugar de eso, calcuaremos *perfiles columna*. Comparamos cada una de las columnas con la
columna marginal (la tabla de tipo de estilo de té):

```{r, message = FALSE}
tabla <- te %>% group_by(how, price) %>% tally() %>%
  group_by(how) %>% 
  mutate(prop_price = (100 * n / sum(n))) %>% 
  group_by(price) %>% 
  mutate(prom_prop = mean(prop_price)) %>% 
  mutate(perfil = (prop_price / prom_prop - 1) %>% round(2))  
```

```{r, echo = FALSE, message = FALSE}
precio_prom <- tabla %>% select(price, prom_prop) %>% unique %>% 
  mutate(promedio = round(prom_prop)) %>% select(price, promedio)
tabla_perfil <- tabla %>%   
  select(how, price, perfil) %>% spread(how, perfil, fill = -1) 
tabla_2 <- tabla_perfil %>% 
  gather(how, prop_price, -price)
if_profile <- function(x){
  any(x < 0) & any(x > 0)
}
marcar <- marcar_tabla_fun(0.25, "red", "black")
tab_out <- tabla_perfil %>% left_join(precio_prom) %>%
  arrange(desc(`tea bag`)) %>% 
  mutate_if(if_profile, marcar) %>% 
  knitr::kable(format = "html", escape = F, digits = 2) %>% 
  kableExtra::kable_styling(bootstrap_options = c( "hover", "condensed"), full_width = FALSE)
tab_out
```

Leemos esta tabla como sigue: por ejemplo, los compradores de té suelto compran té *upscale* (fino) a una
tasa casi el doble (0.98) que el promedio. 

También podemos graficar como:

```{r, fig.width = 6, fig.height = 2, echo = FALSE, message = FALSE}
tabla_ordenada <- tabla %>% ungroup %>% 
  left_join(tabla %>% ungroup %>% filter(how == "tea bag") %>% select(price, perfil_tea = perfil)) %>% 
  mutate(precio = fct_reorder(price, perfil_tea))
g_perfil <- ggplot(tabla_ordenada,
  aes(x = precio, xend = precio, y = perfil, yend = 0, group = how)) + 
  geom_point() + geom_segment() + facet_wrap(~how) +
  geom_hline(yintercept = 0 , colour = "gray")+ coord_flip()
g_props <- ggplot(precio, aes(y = price, x = 1, label = prop, fill = 1)) + 
    scale_fill_gradient(low = "white", high = "white") +
  geom_tile() + geom_text() + theme(axis.text.x = element_blank(), legend.position = "none") +
  xlab(" ") 
g_perfil
```

**Observación**: hay dos maneras de construir la columna promedio: tomando los porcentajes sobre todos los datos,
o promediando los porcentajes de las columnas. Si los grupos de las columnas están desbalanceados, estos promedios
son diferentes.

- Cuando usamos porcentajes sobre la población, perfiles columna y renglón dan el mismo resultado
- Sin embargo, cuando hay un grupo considerablemente más grande que otros, las comparaciones se vuelven vs 
este grupo particular. No siempre queremos hacer esto.




## Comparaciones causales {-}

Los principios del diseño analítico son especialmente importantes
cuando queremos interpretar las comparaciones que hacemos de forma *causal*,
que muchas veces son las que más nos interesan.

Por ejemplo, en el estudio de estereogramas que acabamos de ver, la interpretación que
hacemos final es causal : la descripción verbal reduce el tiempo de fusión en 2/3,
comparado con no recibir la descripción verbal. Esto lo podemos hacer pues 
este caso asignamos al azar el tratamiento de descripción verbal sobre los
individuos del estudio. 

Esta aleatorización elimina diferencias sistemáticas entre los grupos, de forma
que llegar a la interpretación causal es simple: los dos grupos que comparamos son
similares en promedio, las diferencias que puedan tener se deben al azar. No es
necesario (aunque puede ayudar) "refinar" las comparaciones con otras variables
como hicimos en los ejemplos anteriores.

Sin embargo, queda variación no sistemática que también dificulta interpretar
la comparación. En este caso, podemos
usar las herramientas estadísticas de la siguiente sección para evaluar la evidencia
de que los grupos tienen respuestas distintas. Si esa evidencia es fuerte, podemos
adjudicarla a la única diferencia sistemática, que es si recibió la información verbal
o no, y así concluir con una interpretación causal.

- Los experimentos no son la única manera (pero sí la más simple y confiable desde el 
punto de vista analítico) de hacer **comparaciones causales**. 
- Las comparaciones usuales que hicimos en la mayoría de los ejemplos anteriores les podemos llamar **comparaciones predictivas** (más de esto en la sesión 3).

Notése que esto no quiere decir que cuando hacemos comparaciones predictivas usuales debemos dejar de  pensemos en causalidad, sino al contrario:

- Pensamos qué otras variables pueden estar produciendo sistemática o no sistemáticamente las diferencias o contrastes que vemos.
- Pensamos si podemos conseguir o no información de esas variables para refinar nuestras comparaciones.
- Pensamos cómo podemos usar la información que tenemos para hacer comparaciones más interpretables.
- Discutimos por qué pueden ser deficientes nuestras comparaciones.

Lectura adicional: el capítulo 9 de @GelmanHill tiene una introducción accesible.