Modelos Estadísticos. Grado Biotecnología

Análisis exploratorio de una o dos variables

Abstract

En este tema introducimos los procedimientos para el análisis descriptivo de una o dos variables. En el caso de una variable (de tipo numérico o categórico) estamos interesados en el estudio de la distribución del conjunot de datos muestrales, mientras que cuando trabajamos con dos varaibles estamos interesados en la covariación de una variable respecto de la otra. Introduciremos los procedimientos numéricos y gráficos necesarios apra dichos análisis.

Introducción

Esta unidad mostrará los procedimientos de análisis estadístico descriptivo para el estudio de una o dos variables de tipo numérico y/o categórico. Las situaciones que planteamos son:

Una variable de tipo factor.
Una variable de tipo numérico.
Dos variables categóricas.
Dos variables numéricas.
Una variable categórica y una variable numérica.

Este es el primer paso del EDA y sirve al investigador para plantear las primeras preguntas de interés sobre sus datos. Para ilustrar los procedimientos utilizaremos el conjunto de datos storms de la librería nasaweather. Estos datos son un subconjunto de la base de datos de la NASA sobre los huracanes en el Atlántico Norte (NOAA). Los datos incluyen las posiciones y atributos de 198 tormentas tropicales, medidas cada seis horas durante la vida de la tormenta. Las variables registradas son:

name: Nombre de la tormenta
year, month, day: Año, mes y día del informe de la tormenta
hour: Hora del informe (en UTC)
lat, long: Latitud y longitud de la tormenta
pressure: Presión atmosférica en el centro de la tormenta (en milibares)
wind: Máxima velocidad sostenida de la tormenta (en nudos)
type: Clasificación de la tormenta (Tropical Depression, Tropical Storm, or Hurricane)
seasday: día de la temporada de tormentas

Instalamos la librería y la cargamos

library(nasaweather)

Vemos la estructura del banco de datos

storm <- nasaweather::storms # Guardamos los datos en un nuevo objeto
str(storm)

## Classes 'tbl_df', 'tbl' and 'data.frame':    2747 obs. of  11 variables:
##  $ name    : chr  "Allison" "Allison" "Allison" "Allison" ...
##  $ year    : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
##  $ month   : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ day     : int  3 3 3 3 4 4 4 4 5 5 ...
##  $ hour    : int  0 6 12 18 0 6 12 18 0 6 ...
##  $ lat     : num  17.4 18.3 19.3 20.6 22 23.3 24.7 26.2 27.6 28.5 ...
##  $ long    : num  -84.3 -84.9 -85.7 -85.8 -86 -86.3 -86.2 -86.2 -86.1 -85.6 ...
##  $ pressure: int  1005 1004 1003 1001 997 995 987 988 988 990 ...
##  $ wind    : int  30 30 35 40 50 60 65 65 65 60 ...
##  $ type    : chr  "Tropical Depression" "Tropical Depression" "Tropical Storm" "Tropical Storm" ...
##  $ seasday : int  3 3 3 3 4 4 4 4 5 5 ...

Dado que las variables year y month no vienen definidas como factores, el primer paso es convertirlas en factores. En este caso vamos a utilizar una versión diferente de la vista en el tema anterior para convertir varaibles enteras a factores.

# Pirmero el año
# Creamos el factor
storm$year_f <- factor(storm$year)
# Asignamos los niveles
levels(storm$year_f) <- as.character(1995:2000)

# Ahora los meses
# Creasmos el factor
storm$month_f <- factor(storm$month)
# Asignamos los niveles
levels(storm$month_f) <- c("June","July","August","September","October","November","December")

# Veamos como queda el banco de datos
str(storm)

## Classes 'tbl_df', 'tbl' and 'data.frame':    2747 obs. of  13 variables:
##  $ name    : chr  "Allison" "Allison" "Allison" "Allison" ...
##  $ year    : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
##  $ month   : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ day     : int  3 3 3 3 4 4 4 4 5 5 ...
##  $ hour    : int  0 6 12 18 0 6 12 18 0 6 ...
##  $ lat     : num  17.4 18.3 19.3 20.6 22 23.3 24.7 26.2 27.6 28.5 ...
##  $ long    : num  -84.3 -84.9 -85.7 -85.8 -86 -86.3 -86.2 -86.2 -86.1 -85.6 ...
##  $ pressure: int  1005 1004 1003 1001 997 995 987 988 988 990 ...
##  $ wind    : int  30 30 35 40 50 60 65 65 65 60 ...
##  $ type    : chr  "Tropical Depression" "Tropical Depression" "Tropical Storm" "Tropical Storm" ...
##  $ seasday : int  3 3 3 3 4 4 4 4 5 5 ...
##  $ year_f  : Factor w/ 6 levels "1995","1996",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ month_f : Factor w/ 7 levels "June","July",..: 1 1 1 1 1 1 1 1 1 1 ...

Se han creado dos nuevas varaibles que identifican los factores de año y mes.

A lo largo de cada una de las secciones siguientes se irán introduciendo las funciones y procedimientos necesarias para cada análisis.

Antes de comenzar cargamos el resto de librerías de trabajo

library(tidyverse)
library(stringr)
library(forcats)
library(lubridate)
library(magrittr)
library(broom)
library(datasets)

Casi todo el material se ha obtenido de los manuales: Wickham and Grolemund (2016).

Una variable factor

En esta sección se considera cómo explorar la distribución de una variable categórica. Se presentan las descriptivas básicas y visualizaciones que son apropiadas para las variables categóricas. Para ejemplificar los procedimientos utilizaremos la variable type que identifica el tipo de tormenta.

Explorar variables categóricas es generalmente más simple que trabajar con variables numéricas porque tenemos menos respuestas posibles. La primera pregunta que nos debemos plantear es is la variables categórica es de tipo nominal u ordinal. Esto tiene un efecto relevante a la hora de presentar y visualizar la información. La variable type es una variable categórica de tipo ordinal debido a que las tormentas son clasificadas según su virulencia.

Resúmenes numéricos

Cuando calculamos resúmenes de variables categóricas, intentamos describir la distribución muestral de la variable mediante el recuento de ocurrencias de cada una de las posibles respuestas de la variable. Estos recuentos es lo que denominamos en estadística frecuencias absolutas. Necesitamos entender qué categorías son comunes y cuáles son raras. Asociado a dichas frecuencias absolutas podemos obtener las frecuencias relativas, o más frecuentemente los porcentajes de cada una de las posibles categorías de la variable.

Si \(F_i\) es la frecuencia absoluta de la categoría \(i\) de la variable y \(n\) es el tamaño de la muestra, se define la frecuencia relativa de la categoría \(i\) (\(f_i\)) como: \[f_i =\frac{F_i}{n}\] El porcentaje de la categoría \(i\) es simplemente \(100 * f_i\).

La forma más sencilla para obtener la tabla de frecuencias absoluta asociada con una variable categórica es la función table().

table(storm$type)

## 
##       Extratropical           Hurricane Tropical Depression 
##                 412                 896                 513 
##      Tropical Storm 
##                 926

Para obtener los porcentajes en lugar de los conteos podemos utilizar le código siguiente

type_counts <- table(storm$type)
round(type_counts / sum(type_counts),2) # Redondeamos a dos decimales

## 
##       Extratropical           Hurricane Tropical Depression 
##                0.15                0.33                0.19 
##      Tropical Storm 
##                0.34

Sin embargo, hay otra opción un poco más compleja pero que nos servirá de utilidad para el resto de esta unidad y que pasamos a mostrar a continuación. Es una combinación del operador de anidamiento %>%, de la función de agrupación (group_by()), de la función de resumen (summarise()), y de la función de conteos (n()).

# banco de datos
tabla_tipo <- storm %>%
  # agrupamos por la variable factor
  group_by(type) %>%  
  # resumimos contando el número de casos de cada nivel del factor
  summarise(n=n())
# Para calcular los porcentajes
mutate(tabla_tipo,percent=round(100*n/sum(n),2))

Aunque esta forma es más costosa cuando tenemos únicamente una variable nos resultará de más utilidad cuando deseemos realizar análisis que involucren un mayor número de variables.

Podemos ver que mayoritariamente se han producido Tormentas Tropicales (33.71%) y Huracanes (32.62%). Estos dos tipos suman más de dos tercios de los ecentos registrados. Por contra los eventos menos abundantes han sido las Depresiones Tropicales (18.67%) y las Extratropicales (15.00%).

Las funciones más habituales que se utilizan con summarise() son:

Localización: mean() (media), median() (mediana)
Escala: sd() (desviación típica), IQR() (rango intercuartílico)
Rango: min() (mínimo), max() (máximo), quantile() (cuantil)
Posición: first() (primero), last() (último), nth() (posición n-ésima)
Conteo: n() (número de casos), n_distinct() (número de casos distintos)

La mayoría se usan exclusivamente para variables de tipo numérico.

Visualización gráfica

En este apartado veremos como representar los datos de los conteos de una variable categórica mediante la función ’ggplot`. Esta función permite realizar casi cualquier tipo de gráfico que podamos imaginar. En este punto vamos a ir presentando diferentes parámetros de dicha función para ir familiarizándonos con su uso.

Gráfico barras

La herramienta gráfica más común utilizada para resumir una variable categórica es un gráfico de barras. Un gráfico de barras (o gráfico de barras) es un gráfico que presenta resúmenes de datos agrupados con barras rectangulares. La longitud de las barras es proporcional a los valores que representan. Al resumir una sola variable categórica, la longitud de las barras debe mostrar los recuentos brutos o las proporciones de cada categoría.

Para realizar este gráfico con la función ggplotnecesitamos identificar el conjunto de datos sobre el que vamos a trabajar y la variable que queremos representar:

# Configuramos el gráfico identificando los datos y la variable 
# El comando aes() indica la varaible a representar y el eje donde se situarán los valores de dicha varaible
bar_plt <- ggplot(storm, aes(x = type)) 
# Seleccionamos el tipo gráfico
bar_plt <- bar_plt + geom_bar()
# Representamos el gráfico
bar_plt

Esta es la misma información de resumen que producimos usando la función de tabla, solo que ahora se presenta en forma gráfica. Podemos personalizar este gráfico de barras si es necesario con funciones como xlab y ylab, y configurando varias propiedades dentro de geom_bar. Por ejemplo:

# Retocamos las barras para que aparezcan en azul y con un ancho inferior
# En este caso no almacenamos el gráfico sino que lo ejecutamos directamente
ggplot(storm, aes(x = type)) + 
  geom_bar(fill = "blue", width = 0.7) + 
  xlab("Tipo de Tormenta") + ylab("Número de observaciones")

Como podemos ver tanto en las tablas obtenidas como en los dos gráficos precedentes la escala del tipo de tormenta no está ordenada, es decir, no la tenemos graduada por la relevancia de la tormenta. Veamos como podemos hacer esto e integrarlo en el gráfico:

# Creamos un vector con el orden predefinido
ords <- c("Tropical Depression", "Extratropical", "Tropical Storm", "Hurricane")
# Generamos el gráfico indica que el eje x tiene escala dada por el vector ordenado
ggplot(storm, aes(x = type)) + 
  geom_bar(fill = "blue", width = 0.7) + 
  scale_x_discrete(limits = ords) +
  xlab("Tipo de Tormenta") + ylab("Número de observaciones")

Ahora el gráfico si está ordenado con la escala adecuada y resulta más fácil cuantificar la relevancia de las tormentas más importantes.

También podemos intercambiar las filas por las columnas para una mejor visualización de las etiquetas de la variable categórica. Para ello utilizamos el parámetro coord_flip():

ggplot(storm, aes(x = type)) + 
  geom_bar(fill = "blue", width = 0.7) + 
  scale_x_discrete(limits = ords) +
  coord_flip() + 
  xlab("Tipo de Tormenta") + ylab("Número de observaciones")

Por último utilizamos la función theme_bw() para configurar un fondo blanco para el gráfico

ggplot(storm, aes(x = type)) + 
  geom_bar(fill = "blue", width = 0.7) + 
  scale_x_discrete(limits = ords) +
  coord_flip() + 
  xlab("Tipo de Tormenta") + ylab("Número de observaciones")+
  theme_bw()

Otras posibilidades para los temas son theme_classic(), theme_dark(), theme_grey(), theme_light().

Repreesntamos ahora los porcentajes asociados a cada categoría en lugar de los conteos:

ggplot(storm, aes(x = type)) + 
  geom_bar(aes(y = ..prop.. , group = 1),fill = "blue", width = 0.7) + 
  scale_y_continuous(labels = scales::percent) +
  coord_flip() + 
  xlab("Tipo de Tormenta") + ylab("Porcentaje")+
  theme_bw()

Una variable numérica

En esta sección se considera cómo explorar la distribución de una variable numérica Se presentan las descriptivas básicas y visualizaciones que son apropiadas para las variables de este tipo. Para ejemplificar los procedimientos utilizaremos las variables wind y pressure.

Resúmenes numéricos

Hasta ahora hemos estado describiendo las propiedades de las distribuciones de muestra en términos muy generales, usando frases como “valores más comunes” y “el rango de los datos” sin decir realmente lo que queremos decir. Los estadísticos han ideado términos específicos para describir este tipo de propiedades, así como diferentes estadísticas descriptivas para cuantificarlas. Los dos que más importan son la tendencia central y la dispersión:

Una medida de tendencia central describe un valor típico (‘central’) de una distribución de datos. La medida de localización más extendida es la media aritmética de una muestra. Hay muchas medidas diferentes de tendencia central, cada una con sus propios pros y contras. Entre estos, la mediana es la que se usa con mayor frecuencia en los análisis exploratorios ya que es le valor que nos divide la muestra de dos partes iguales situando el 50% de los datos a cada lado de ese valor.
Una medida de dispersión describe cómo se distribuye una distribución. Las medidas de dispersión cuantifican la variabilidad o dispersión de una variable con respecto al promedio de los datos. Si una distribución está más dispersa que otra, significa que, en cierto sentido, abarca una gama más amplia de valores. Lo que esto significa en la práctica depende del tipo de medida con la que estamos trabajando. Las medidas de dispersión más habituales son la varianza y su raíz cuadrada, la desviación estándar.

Medidas de tendencia central

Hay dos estadísticos que se utilizan generalmente para describir la tendencia central de la distribución de los datos muestrales de un variable numérica. De ahora en adelante denotamos por \(n\) al tamaño muestral y \(x_1, x_2,...,x_n\) los valores muestrales de la variable que deseamos estudiar.

La media muestral es la medida de tendencia muestral por excelencia. La definición matemática de la media muestral viene dada por:

\[\bar{x} = \frac{\sum_{i=1}^n x_i}{n}\] Para obtener la media utilizamos la función mean()

mean(storm$wind)

## [1] 54.68329

mean(storm$pressure)

## [1] 989.8238

Esto nos dice que la media de la velocidad del viento es de 55 mph y que la media de la presión es de 989.82 milibares. ¿Como podemos interpretar esos resultados?

Una limitación de la media aritmética es que se ve afectada por la forma de la distribución de los datos. Es muy sensible a los extremos de una muestra. Esta es la razón por la cual, por ejemplo, no tiene mucho sentido mirar el ingreso medio de los trabajadores en un país para tener una idea de lo que gana una persona “típica”. La distribución del ingreso es muy asimétrica, y los pocos que tienen la suerte de ganar salarios muy buenos tienden a cambiar la media hacia arriba y superar cualquier cosa que sea realmente “típica”. La media de la muestra también se ve fuertemente afectada por la presencia de ‘valores atípicos’ o valores extremos, es decir, valores inusualmente grandes o pequeños en una muestra.

Debido a que la media muestral es sensible a la forma de una distribución y la presencia de valores atípicos, a menudo se prefiere una segunda medida de tendencia central: la mediana de la muestra. La mediana de una muestra es el número que separa los datos en dos subgrupos (la mitad superior de la mitad inferior). Podemos calcular la mediana muestral en R con la función median():

median(storm$wind)

## [1] 50

median(storm$pressure)

## [1] 995

Estos resultados indican que el 50% de de los registros muestran una valor de viento inferior a 50 mph. De la misma forma el 50% de los datos muestran una valor de la presión inferior a 995 milibares.

Otras medidas de localización son el mínimo, el máximo y los percentiles. Los percentiles son los valores que dividen en la muestra según el valor del percentil solicitado. Si solicitamos el percentil 20, se separa la muestra en dos subconjuntos dejando el 20% en un grupo y e 80% en el otro. Los percentiles más habituales son lo denominados primer y tercer cuartil que corresponden a los percentiles 25 y 75 respectivamente.

Para obtener todas las medidas de localización de una vez utilizamos la función summary()

summary(storm$wind)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.00   35.00   50.00   54.68   70.00  155.00

summary(storm$pressure)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   905.0   980.0   995.0   989.8  1004.0  1019.0

Se pude ver que el 25% de las observaciones muestran una valor de viento inferior a 35mph y un 25% un valor de viento superior a 70mph. Con respecto a la presión atmosférica tenemos que el 25% de las observaciones tienen un valor inferior a 980 milibares, y un 25% tienen un valor superior a 1004 milibares.

Medidas de dispersión

Hay muchas maneras de cuantificar la dispersión de un conjunto de datos muestrales de una variable numérica. Los valores más importantes desde el punto de vista estadístico son la varianza muestral y la desviación estándar. La varianza muestral \(s^2\) es “la suma de las desviaciones cuadradas” (es decir, las diferencias) de cada observación con respecto a la media de la muestra, dividida por el tamaño de la muestra menos uno. La desviación típica es la raíz cuadrada de la varianza muestral. Las definiciones matemáticas de ambas cantidades son: \[s^2 = \frac{\sum_{i=1}^n (x_i-\bar{x})^2}{n-1}\] \[s = \sqrt{s^2}\] Las funciones de R para calcular ambas cantidades son var() para la varianza y sd() para la desviación típica.

var(storm$wind);sd(storms$wind)

## [1] 668.1444

## [1] 26.21387

var(storm$pressure);sd(storms$pressure)

## [1] 349.4912

## [1] 19.51678

¿Qué significa ese número en realidad? Las variaciones son siempre no negativas. Una pequeña varianza indica que las observaciones tienden a ser cercanas a la media (y una a la otra), mientras que una alta varianza indica que las observaciones están muy dispersas. Una varianza de cero solo ocurre si todos los valores son idénticos. Sin embargo, es difícil interpretar si una varianza muestral es realmente “pequeña” o “grande” porque el cálculo involucra desviaciones al cuadrado. Por ejemplo, cambiar la escala de medición de una variable por 10 implica un cambio de 100 veces (102) en la varianza.

La varianza es una cantidad importante en las estadísticas que aparece una y otra vez. Muchas herramientas estadísticas comunes usan cambios en la varianza para comparar formalmente qué tan bien diferentes modelos describen un conjunto de datos. Sin embargo, es muy difícil interpretar las variaciones, por lo que rara vez las utilizamos en el trabajo exploratorio. Para expresar la variabilidad en la misma escala de la variable original utilizamos la desviación típica. En este caso se puede observar una mayor variabilidad en la variable presión atmosférica (desviación típica de 25.84 milibares) que en la variable de viento (18.69 mph).

La desviación estándar de la muestra no está exenta de problemas. Al igual que la media muestral, es sensible a la forma de la distribución de los datos y a la presencia de valores atípicos. Una medida de dispersión más robusta para este tipo de características es el rango intercuartílico, definida como la diferencia entre el percentil 75 (tercer cuartil) y el percentil 25 (primer cuartil): \[IQR = Q_3 - Q_1\] Obviamente, cuanto más dispersos estén los datos, mayor será el IQR. La razón por la que preferimos usar IQR para medir la dispersión es que solo depende de los datos en el “medio” de una distribución de muestra. Esto lo hace robusto a la presencia de valores atípicos. Podemos usar la función IQR() para calcular el rango intercuartílico:

IQR(storm$wind)

## [1] 35

IQR(storm$pressure)

## [1] 24

Visualización gráfica

El gráfico por excelencia para una variable de tipo numérico es el denominado histograma. El histograma es una representación mediante barras de la distribución de los datos. Para construir las barras se divide el rango de la variable numérica en un conjunto fijo de intervalos disjuntos y se contabiliza el número de datos que quedan dentro de cada uno de ellos (altura del gráfico de barras). Es un gráfico muy interesante porque representa de forma bastante precisa, ajustando el número de intervalos, la distribución del conjunto de datos pudiéndose observar su dispersión y/o asimetria.

Para realizar este gráfico utilizamos parámetros específicos dentro de la función ggplot()

ggplot(storm, aes(x = pressure)) +
   geom_histogram()

Podemos ver que la mayoría de los datos se sitúan por encima de 915 milibares, hay una clara asimetría hacia los valores más grandes y existe una gran dispersión entre en el conjunto de valores observados.

Veamos como introducir diferentes parámetros en el gráfico anterior. el más importante es el que hace referencia a los intervalos asociado con el histograma (parámetro binwith)

ggplot(storm, aes(x = pressure)) + 
  geom_histogram(binwidth = 8, fill = "steelblue") + 
  xlab("Presión Atmosférica (mbar)") + ylab("Conteo")

Veamos ahora el histograma de la variable wind

ggplot(storm, aes(x = wind)) + 
  geom_histogram(binwidth = 8, fill = "steelblue") + 
  xlab("Velocidad del viento (mph)") + ylab("Conteo")

El otro gráfico habitual para una variable de tipo numérico es el llamado gráfico de cajas. En este gráfico se representa mediante una caja la media (linea central de la caja), el percentil 75 (línea superior de la caja), y el percentil 25 (linea inferior de la caja). También se representan el valor máximo (percentil 75 + 1.5IQR) y mínimo (percentil 25 - 1.5IQR), así como los caracterizados como valores extremos (punto fuera de la caja). Veamos como realizar este gráfico mediante el parámetro geom_boxplot().

ggplot(storm, aes(x = factor(1),y = pressure)) +
   geom_boxplot()

Añadimos parámetros (título de ejes, color caja, selección de valores extremos, y fondo blanco)

ggplot(storm, aes(x = factor(1),y = pressure)) +
   geom_boxplot(fill = "orange",outlier.colour = "red", outlier.shape = 1) +
   scale_x_discrete(name = " ") + scale_y_continuous(name = "Presión Atmosférica (mbar)") +
   theme_bw()

Podemos ver que la media se sitúa muy próxima a los 1000 mbar y la gran cantidad de valores extremos en la parte baja de la distribución. Además la caja es muy estrecha indicando la poca variabilidad en los datos, lo que viene corroborado también por la proximidad de los percentiles 25 y 75.

Ahora le gráfico para la variable wind

ggplot(storm, aes(x = factor(1),y = wind)) +
   geom_boxplot(fill = "orange",outlier.colour = "red", outlier.shape = 1) +
   scale_x_discrete(name = " ") + scale_y_continuous(name = "Velocidad del Viento (mph)") +
   theme_bw()

¿Cómo interpretamos este gráfico?

Dos variables categóricas

Explorar numéricamente las asociaciones entre pares de variables categóricas no es tan simple como el caso de una variable. La pregunta general que debemos abordar es, “¿las diferentes combinaciones de categorías parecen estar sub o sobre representadas?”. Necesitamos entender qué combinaciones son comunes y cuáles son raras. Lo más simple que podemos hacer es construir una tabulación cruzada del número de ocurrencias de cada combinación de niveles de ambas variables. La tabla resultante se llama tabla de contingencia.

En cuanto a la representación gráfica la opción más habitual pasa por representar los conteos mediante gráficos de barras que representan de forma conjunta la información de ambas variables.

Para ejemplificar los cálculos y gráficos utilizaremos las variables month y type.

Resúmenes numéricos

Para obtener tanto la tabla de frecuencias como los porcentajes asociados vamos a utilizar la función de anidamiento que hemos introducido anteriormente, de forma similar a los cálculos realizados para una variable categórica.

Calculamos la tabla de frecuencias y porcentajes

# Realizamos los cálculos ahora
tabla_dbl <- storm %>% 
  count(month_f,type) %>% 
  mutate(porcentaje=round(100*n/sum(n),2))
tabla_dbl

Para poder visualizar los resultados en formato de tabla haremos uso de la función spread(). En primer lugar obtenemos la tabla de frecuencias de ambas variables:

tabla_conteos <- select(tabla_dbl,month_f,type,n)
spread(tabla_conteos, key = type, value = n)

Lo más relevante en estas tablas es buscar los valores más altos y más bajos para detectar las combinaciones más y menos habituales. En este caso la combinación más habitual son los huracanes en el mes de septiembre, seguida de los huracanes en agosto. ¿Cuáles son las combinaciones menos relevantes?

Para apreciar mejor este tipo de relaciones representamos la tabla de porcentajes:

tabla_porcentajes <- select(tabla_dbl,month_f,type,porcentaje)
spread(tabla_porcentajes, key = type, value = porcentaje)

Podemos ver que la combinación Huracán en el mes de septiembre representa el 13.94% de los datos observados, mientras que la segunda combinación más frecuente (huracanes - agosto) representa el 10.92%. Juntos representan casi una cuarta parte del total de datos observados. La combinación menos frecuente (0.04%) corresponde con las tormentas tropicales en el mes de diciembre.

Visualización gráfica

Los gráficos de barras se pueden usar para resumir la relación entre dos variables categóricas. La idea básica es producir una barra separada para cada combinación de categorías en las dos variables. La longitud de estas barras es proporcional a los valores que representan, que son los recuentos brutos o las proporciones en cada combinación de categorías. Esta es la misma información que se muestra en una tabla de contingencia. El uso de ggplot2 para mostrar esta información no es muy diferente de producir un gráfico de barras para resumir una única variable categórica.

Tomamos las variables type y year para mostrar su funcionamiento. Ordenamos la variable type para mostrar los gráficos por orden de importancia de la tormenta.

# Creamos un vector con el orden predefinido
ords <- c("Tropical Depression", "Extratropical", "Tropical Storm", "Hurricane")
# Generamos el gráfico indica que el eje x tiene escala dada por el vector ordenado
ggplot(storm, aes(x = type, fill = year_f)) + 
  geom_bar() + 
  scale_x_discrete(limits = ords) +
  xlab("Tipo de Tormenta") + ylab("Número de observaciones")

Esto se llama un gráfico de barras apiladas. Cada tipo de tormenta tiene su propia barra, y cada barra se ha dividido en diferentes segmentos de colores, cuya longitud está determinada por el número de observaciones asociadas con cada año. Podemos apreciar si para un mismo tipo de tormenta la ocurrencia de un tipo de tormenta es similar o no. Por ejemplo, para los huracanes podemos apreciar que los años 1996, 1998, 1999, y 2000 tienen un número similar de ocurrencias.

Un problema con este tipo de gráfico es que puede ser difícil detectar asociaciones entre las dos variables categóricas. Si queremos saber cómo están asociados, a menudo es mejor trazar los recuentos para cada combinación de categorías una al lado de la otra. Utilizamos la opción dodge en la función geom_bar:

ggplot(storm, aes(x = type, fill = year_f)) + 
  geom_bar(position = "dodge") + 
  scale_x_discrete(limits = ords) +
  labs(x = "Tipo de tormenta", y = "Número de observaciones", fill = "Año")

¿Qué ventajas o desventajas aprecias en este gráfico frente al de barras apiladas?

Veamos el gráfico de porcentajes:

ggplot(storm, aes(x = type, fill = year_f)) + 
  geom_bar(aes(y = ..prop.. , group = year_f),position = "dodge") + 
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Tipo de tormenta", y = "Porcentaje", fill = "Año")

Otra opción es realizar un mapa de intensidad de cada una de las combinaciones e ambos factores. Utilizamos la función geom_tile().

storm %>%
    count(type,year_f) %>%
    ggplot(mapping = aes(x = type, y = year_f)) + 
        geom_tile(mapping = aes(fill = n)) + 
        scale_x_discrete(limits = ords) +
        labs(x = "Tipo de tormenta", y = "Año", fill = "n")

En este caso las casillas en tonos más claros corresponden con las combinaciones de niveles con una mayor ocurrencia. Las más abundantes se concentran en 1995 y los tipos de tormenta más graves. Por otro lado el año 1997 es el que ha registrado menor número de tormentas.

Categórica vs Numérica

El objetivo en este tipo de situaciones es comparar cada una de las distribución de la variable numéricas que surgen al segmentar los valores por cada uno de los niveles de la variable categórica. Creamos tanto grupos como niveles del factor.

Resúmenes numéricos

Los resúmenes numéricos se pueden construir tomando las diversas ideas que hemos explorado para las variables numéricas (medias, medianas, etc.) y aplicándolas a subconjuntos de datos definidos por los valores de la variable categórica. Utilizamos la función group_by para agrupar los valores de la variable numérica y realizar todos los cálculos. Utilizamos la variable type como categórica y las variables wind y pressure como numéricas.

storm %>% 
  group_by(type) %>% 
  summarise(
    conteo = n(),
    m_wind = round(mean(wind, na.rm = TRUE),2),
    dt_wind = round(sd(wind, na.rm = TRUE),2),
    m_pressure = round(mean(pressure, na.rm = TRUE),2),
    dt_pressure = round(sd(pressure, na.rm = TRUE),2)
  )

La opción `na.rm = TRUE indica que no se deben tener en cuenta los valores perdidos para el cálculo solicitado. Se puede ver como el tipo Huracanes tiene la media más alta de velocidad del viento y la más baja de presión atmosférica. Además en ambos casos la variabilidad observada es la mayor de todos los tipos de tormenta. La depresión tropical es justo el caso contrario. Evidentemente los datos observados corresponden con la clasificación de tormenta establecida desde el inicio. Podemos obtener cualquier otra medida de tipo descriptivo de las que vimos en el apartado correspondiente a una variable numérica, ya que al dividir por los niveles el factor es como si tuviéramos cuatro nuevas variables.

Otra forma de realizar los cálculos anteriores escribiendo menos código es utilizando la función summarise_each que permite calcular un conjunto de funciones sobre un conjunto de variables.

storm %>% 
  group_by(type) %>% 
  summarise_each(
    funs(mean,sd),wind,pressure)

Visualización gráfica

Para la visualización gráfica de las relaciones entre una variable categórica y una numérica tenemos diferentes opciones: a) gráfico de densidad gráfico de cajas, y gráficos comparativos matriciales.

El gráfico de densidad representa la distribución del conjunto de datos muestrales. Con este gráfico se puede apreciar claramente el rango de valores y su concentración (altura de la curva de densidad). Para su obtención sólo debemos fijar el parámetro de suavizado, bw, que nos indica el grado de información que debemos utilizar para obtener dicha densidad. Valores pequeños dan curvas poco suavizadas y valores grandes dan curva suavizadas. Siempre es necesario un pequeño ajuste para obtener el valor más adecuado. Con este gráfico resulta muy fácil comparar el comportamiento de los diferentes grupos como veremos en los gráficos siguientes.

Comenzamos con la variable velocidad del viento:

ggplot(storm, aes(x = wind)) + 
  geom_density(aes(colour = type),bw = 3, na.rm = TRUE) + 
  labs(x = "Velocidad del viento (en mph)", y = "Densidad")

El resultado son cuatro curvas de densidad (una por cada top de tormenta) donde los más destacable es que el rango de valores posibles para la velocidad del viento es sólo diferente para el tipo huracanes, ya que su curva de densidad se encuentra desplazada respecto del resto de densidades. También se puede ver que las depresiones tropicales tienen una menor dispersión (menor rango de valores) los que provoca que su curva sea más puntiaguda. Cuanto mayor sea la dispersión mas amplia sera la densidad obtenida. Si la densidad es simétrica la media se sitúa en el punto medio, mientras que el valor asociado con el punto más alto de la densidad es lo que denominamos moda.

Ahora con la variable presión atmosférica:

ggplot(storm, aes(x = pressure)) + 
  geom_density(aes(colour = type),bw = 3, na.rm = TRUE) + 
  labs(x = "Velocidad del viento (en mph)", y = "Densidad")

El comportamiento de la variable presión atmosférica es muy similar al de la velocidad del viento pero en la parte izquierda del rango de valores. Los huracanes tienen la menor presión atmosférica y se distinguen del resto de tipos de tormenta. También muestra una mayor variabilidad que el resto de tipos.

La visualización más común para explorar las relaciones entre una variable categórica y otra numérica es el diagrama de cajas. Cada diagrama de caja consiste en:

Una casilla que se extiende desde el percentil 25 de la distribución hasta el percentil 75, una distancia conocida como rango intercuartílico (IQR). En el medio del recuadro hay una línea que muestra la mediana, es decir, el percentil 50 de la distribución. Estas tres líneas le dan una idea de la extensión de la distribución y si la distribución es simétrica o no respecto a la mediana o sesgada hacia un lado.
Puntos visuales que muestran observaciones que caen más de 1,5 veces el IQR desde cualquier borde de la caja. Estos puntos remotos son inusuales, por lo que se representan de forma individual.
Una línea (o bigote) que se extiende desde cada extremo de la caja y va al el punto más lejano no atípico en la distribución.

Veamos los ejemplos:

ggplot(storm, aes(x = type, y = wind)) + 
  geom_boxplot() + 
  scale_x_discrete(limits = ords) +
  labs(x = "Tipo de Tormenta", y = "Velocidad del viento (en mph)")

En este tipo de gráficos se esta interesado en dos aspecto fundamentales: * Estudiar la variabilidad dentro de cada grupo viendo la altura de la caja (IQR). * Comparar el comportamiento de cada grupo observando si las cajas quedan a alturas superpuestas.

En este caso podemos ver que la variabilidad más grande se produce en el tipo huracanes y la más pequeña en las depresiones tropicales. Se aprecian diferencias entre las variabilidades de los grupos. Por otro lado, las cajas correspondientes a las tormentas extra tropicales y tormentas tropicales quedan a una misma altura indicando que tienen valores similares para la velocidad del viento. Este gráfico nos da unas primeras indicaciones claras para los procedimientos de comparaciones de medias que estudiaremos más adelante.

ggplot(storm, aes(x = type, y = pressure)) + 
  geom_boxplot() + 
  scale_x_discrete(limits = ords) +
  labs(x = "Tipo de Tormenta", y = "Presión atmosférica (en milibares)")

En este caso observamos que los tipos extra tropical y Huracanes muestran variabilidades similares, mientras que los otros dos tipos también muestran variabilidades similares. En cuanto a la comparación de los grupos se observa que el único con un comportamiento diferente son los huracanes. Los otros tres tipos muestran cajas que se podrías solapar mostrando una mayor igualdad en los valores de presión atmosférica.

Los gráficos matriciales o de facetas permiten representar mediante múltiples gráficos la información de una variable numéricas con respecto a los niveles de la variable categórica. Se pretende de esta forma comprobar la forma de la distribución de los datos de forma similar al gráfico de densidad pero realzando un gráfico para cada nivel. Aunque resultan más útiles cuando trabajamos con más de dos variables, se introducen aquí para ir conociendo su funcionamiento en ejemplos sencillos.

Comenzamos realizando un histograma independiente. Para que resulte más fácil visualizar el gráfico introducimos un orden asociado con los valores de la variable numérica, comenzando con el nivel que tiene valores en esa variable más pequeños, y finalizando con la que tiene los valores más grandes.

# Creamos un nuevo factor ordenado de acurdo a la variable que estamos midiendo
storm$type2 <- reorder(storm$type, storm$wind)
# Creamos el gráfico
ggplot(storm, aes(x = wind))  +
  geom_histogram(binwidth = 5) + 
  xlab("Velocidad del viento (mph)") +
  ylab("Conteo") +
  facet_wrap(~ type2, ncol = 1)

También podríamos realizar el gráfico de densidad

# Creamos un nuevo factor ordenado de acurdo a la variable que estamos midiendo
storm$type2 <- reorder(storm$type, storm$wind)
# Creamos el gráfico
ggplot(storm, aes(x = wind))  +
  geom_density(bw = 3) + 
  xlab("Velocidad del viento (mph)") +
  ylab("Densidad") +
  facet_wrap(~ type2, ncol = 1)

En ambos casos las interpretaciones son similares a las que se hicieron con el gráfico conjunto de densidad.

Podemos cambiar la configuración cambiando el número de columnas:

# Creamos un nuevo factor ordenado de acurdo a la variable que estamos midiendo
storm$type2 <- reorder(storm$type, storm$wind)
# Creamos el gráfico
ggplot(storm, aes(x = wind))  +
  geom_density(bw = 3) + 
  xlab("Velocidad del viento (mph)") +
  ylab("Densidad") +
  facet_wrap(~ type2, ncol = 2)

Aunque el gráfico se visualiza mejor también resulta más difícil la comparación entre todos los niveles.

Si deseamos cambiar la situación de las etiquetas del factor podemos utilizar la opción facet_grid.

# Creamos un nuevo factor ordenado de acurdo a la variable que estamos midiendo
storm$type2 <- reorder(storm$type, storm$wind)
# Creamos el gráfico
ggplot(storm, aes(x = wind))  +
  geom_density(bw = 3) + 
  xlab("Velocidad del viento (mph)") +
  ylab("Densidad") +
  facet_grid(type2 ~ .)

Esta opción nos resultará de mayor utilidad cuando tengamos que representar dos factores ya que se podrá situar en la filas uno de los factores y el otro en las columnas con la opción facet_grid(factor1 ~ factor2).

Dos variables numéricas

Los estadísticos han ideado varias formas diferentes de cuantificar la asociación entre dos variables numéricas en un banco de datos. Las medidas más comunes medidas buscan calcular algún tipo de coeficiente de asociación Los términos “asociación” y “correlación” están estrechamente relacionados; tanto que a menudo se usan indistintamente. La más habitual es la correlación lineal que cuantifica el grado de asociación lineal entre dos variables de tipo numérico. Para ejemplificar nuestros cálculos y gráficos utilizaremos las variables wind y preassure.

Resúmenes numéricos

La medida de correlación más utilizada es el coeficiente de correlación lineal de Pearson. El coeficiente de correlación de Pearson cuantifica el grado de asociación entre las variables en la escala -1 a 1, donde -1 indica una relación inversa (cuando una variable crece la otra decrece) y 1 indica una relación directa (cuando una crece la otra también lo hace. Valores próximo a cero indican que no hay asociación lineal entre las variables analizadas. Este coeficiente es la base para plantear lo que denominaremos más adelante los modelos de regresión lineal simple.

La definición formal del coeficiente de correlación de Pearson (\(\rho\)) viene dada por: \[\rho = \frac{1}{n-1}\sum_{i=1}^n \frac{(x_i - \bar{x})(y_i - \bar{y})}{s_x s_y}\] donde \(x_i\), \(y_i\) son las observaciones de la variable \(x\) e \(y\) respectivamente, \(\bar{x}\), \(\bar{y}\) son las medias muestrales de cada variable, y \(s_x\), \(s_y\) son las desviaciones típica de cada variable. El coeficiente trata de valorar la “covaraición” entre ambas variables, es decir, como afectan los cambios de valores en una variable en los valores de la otra, teniendo en cuenta la propia variabilidad de cada una de ellas.

Para obtener el coeficiente de correlación utilizamos la función cor().

cor(storm$wind,storm$pressure)

## [1] -0.9254911

El coeficiente de correlación resulta negativo, lo que indica que la velocidad del viento tiende a disminuir al aumentar la presión. Al estar próximo a -1 se puede intuir que dicha asociación es muy fuerte. Sin embargo, el coeficiente de correlación de Pearson debe interpretarse con cuidado debido a que está diseñado para medir una relación de tipo lineal, lo que implica que dicho coeficiente será engañoso cuando esta relación sea curva, o incluso peor, en forma de joroba.

¿Qué deberíamos hacer si pensamos que la relación entre dos variables no es lineal? No deberíamos usar el coeficiente de correlación de Pearson para medir la asociación en este caso. En cambio, podemos calcular lo que denominamos correlación de rango. La idea es muy simple. En lugar de trabajar con los valores reales de cada variable, los “clasificamos”, es decir, ordenamos cada variable de menor a mayor y asignamos las etiquetas “primero”, “segundo”, “tercero”, etc. a diferentes observaciones. Las medidas de correlación de rangos se basan en una comparación de los rangos resultantes. Los dos más populares son Spearman’s y Kendall’s. Ambos coeficientes se comportan de una manera muy similar al coeficiente de correlación de Pearson. Toman un valor de 0 si los rangos no están correlacionados, y un valor de +1 o -1 si están perfectamente relacionados.

Podemos calcular ambos coeficientes de correlación de rangos en R usando nuevamente la función cor. Esta vez necesitamos establecer el argumento del método en el valor apropiado: method = "kendall" o method = "spearman".

cor(storm$wind,storm$pressure,method = "kendall")

## [1] -0.7627645

cor(storm$wind,storm$pressure,method = "spearman")

## [1] -0.9025831

Los resultados obtenidos son compatibles con el del coeficiente de Pearson.

Visualización gráfica

Los coeficientes de correlación nos dan una forma simple de resumir las asociaciones entre variables numéricas. Sin embargo, son limitados, porque un solo número nunca puede resumir todos los aspectos de la relación entre dos variables. Es por eso que siempre visualizamos la relación entre dos variables. El gráfico estándar para mostrar asociaciones entre variables numéricas es un diagrama de dispersión, usando ejes horizontales y verticales para trazar dos variables como una serie de puntos. Para realizar este gráfico usamos la opción geom_point()

ggplot(storm, aes(x = wind, y = pressure)) + 
  geom_point() + 
  labs(x = "Velocidad del viento (en mph)", y = "Presión atmosférica (en milibares)")

En el gráfico se puede apreciar la relación de tipo lineal en orden o pendiente decreciente (cuando aumenta el viento disminuye la presión).

El problema de este gráfico es que no podemos apreciar todos los puntos, ya que si tenemos dos observaciones con los mismo valores en ambas variables, estos quedarían superpuestos. Para solucionar esta deficiencia podemos optar por otra versión del gráfico de dispersión que nos permita contabilizar el número de repeticiones cuando estas existan. Este gráfico se obtiene con la opción geom_hex(). Para poder realizarlo es necesario instalar la librería hexbin:

ggplot(storm, aes(x = wind, y = pressure)) + 
  geom_hex(bins = 25) + 
        labs(x = "Velocidad del viento (en mph)", y = "Presión atmosférica (en milibares)", fill = "n")

El parámetro bins segmenta el rango de cada variable en intervalos disjuntos. Lo que se representa es una gráfico de dispersión por intervalos, de forma que cada casilla representa todos los valores que quedan dentro del intervalo conjunto que obtenemos con ambas variables. Se observa que la tendencia se mantiene pero resulta posible ver que valores muestran una mayor o menor ocurrencia. La combinación de valores bajos de viento (< 40 mph) con altos de presión (> 990 mb) son los que más aparecen en el banco de datos.

Otra opción es agrupar una variable continua para que actúe como una variable categórica. Luego se puede usar un gráfico combinado de cajas para representar ambas variables Veamos un ejemplo:

ggplot(data = storm, aes(x = wind, y = pressure)) + 
  geom_boxplot(mapping = aes(group = cut_width(wind, 10))) + 
        labs(x = "Velocidad del viento (en mph)", y = "Presión atmosférica (en milibares)", fill = "n")

La interpretación es similar a la que se realizaba cuando trabajamos con una variable factor y otra numérica. Lo que resulta interesante es que podemos observar los intervalos con un mayor volumen de valores extremos o anómalos (valores de viento entre 30 y 90 mph).

Bibliografía

Wickham, Hadley, and Garrett Grolemund. 2016. R for Data Science. Http://r4ds.had.co.nz/. O’Reilly.