Modelos Estadísticos. Grado Biotecnología

Introducción

En esta unidad se amplían los procedimientos de análisis descriptivos vistos en el tema anterior para estudiar una o dos variables de tipo numérico o categórico al caso de más de dos variables de este tipo. No se hace un barrido a cualquier situación que pueda aparecer sino que se pretende mostrar los casos más habituales. Dichos casos son:

Tres variables categóricas.
Dos variables categóricas y una variable numérica.
Una variable categórica y dos variables numéricas.
Dos variables categóricas y dos variables numéricas.
Tres variables categóricas y dos variables numéricas.

De nuevo utilizaremos el conjunto de datos storms de la librería nasaweather. Recordemos que 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

A lo largo de cada una de las secciones siguientes se irán ampliando (siguiendo los visto en el tema anterior) las funciones y procedimientos necesarios para cada análisis.

Antes de comenzar cargamos las librerías de trabajo

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

Hacemos la copia del banco de datos:

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

Convertimos en factores las variables año y mes.

# Prirmero 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 ...

Tres factores

Los procedimientos numéricos se restringen en este caso a la obtención de la tabla de frecuencias conjunta de las tres variables, mientras que los gráficos se basan en gráficos matriciales donde se consideran gráficos de barras.

A modo de ejemplo vamos a realizar el análisis conjunto de las variables year_f, month_f y type. Calculamos la tabla de frecuencias y porcentajes

# Realizamos los cálculos ahora
tabla_dbl <- storm %>% 
  count(year_f,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().

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

¿Qué conclusiones podemos extraer de estos resultados?

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

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

Como se puede ver en las tablas anteriores al ir aumentando el número de variables resulta más difícil extraer conclusiones a partir de las tablas, y resulta más sencillo interpretar los resultados que aparecen en los gráficos correspondientes.

Para realizar le gráfico combinado de las tres variables categóricas utilizamos un gráfico matricial con dos factores y representamos dentro de cada combinación el gráfico de barras de la otra variable.

ords <- c("Tropical Depression", "Extratropical", "Tropical Storm", "Hurricane")
ggplot(storm, aes(x = type))  +
  geom_bar() + 
  scale_x_discrete(limits = ords) +
  xlab("Tipo de tormenta") +
  ylab("Conteo") +
  facet_grid(year_f ~ month_f)+
  theme(axis.text.x = element_text(angle = 90)) # Para cambiar la orientación de las etiquetas

El gráfico, aún proporcionando la misma información que la tabla de frecuencias, resulta mucho más revelador, ya que se aprecian de forma directa las combinaciones de año - mes en al que no hay datos, y en aquellas donde si los hay se puede ver claramente cual es el tipo de tormenta más predominante.

¿Qué conclusiones podemos extraer a la vista del gráfico?

Dado que la mayoría de los datos se producen entre los meses de agosto y octubre vamos a filtrar los datos para estudiar esas combinaciones únicamente.

storm_meses <- storm %>%
  filter(month_f == c("August","September","October"))
ords <- c("Tropical Depression", "Extratropical", "Tropical Storm", "Hurricane")
ggplot(storm_meses, aes(x = type))  +
  geom_bar() + 
  scale_x_discrete(limits = ords) +
  xlab("Tipo de tormenta") +
  ylab("Conteo") +
  facet_grid(year_f ~ month_f) +
  theme(axis.text.x = element_text(angle = 90)) # Para cambiar la orientación de las etiquetas

Ahora se pueden estudiar con más detalle los meses que concentran un mayor número de tormentas.

Dos factores, Una numérica

En este caso generalizamos el cálculo de medidas de localización y dispersión a esta situación, y analizamos los diferentes gráficos que podemos realizar en esta situación. Utilizamos las variables year_f, type, y wind.

En primer lugar calculamos la tabla de medidas descriptivas:

resumen <- storm %>% 
  group_by(year_f,type) %>% 
  summarise_each(
    funs(mean,sd),wind)
resumen

La tabla es difícil de manejar y por eso creamos la tabla de medias y de desviaciones típicas por separado.

# tabla de medias
resumen_medias <- select(resumen,year_f,type,wind_mean)
spread(resumen_medias, key = type, value = wind_mean)

# tabla de desviaciones típicas
resumen_dt <- select(resumen,year_f,type,wind_sd)
spread(resumen_dt, key = type, value = wind_sd)

Como antes el análisis de estas tablas es más complejo que tratar de representar los datos de forma que se puedan extraer conclusiones de forma más efectiva. Empezamos con el gráfico matricial mezclado con los gráficos 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, color = type2))  +
  geom_density(bw = 3) + 
  xlab("Velocidad del viento (mph)") +
  ylab("Densidad") +
  facet_grid(year_f ~ .)

Ahora realizamos el gráfico de cajas con una estructura similar

# 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 = type2, y = wind))  +
  geom_boxplot() + 
  xlab("Tipo de tormenta") +
  ylab("Velocidad del viento (mph)") +
  facet_grid(. ~ year_f) +
  theme(axis.text.x = element_text(angle = 90))

¿Qué conclusiones podemos extraer de este gráfico?

Otra versión de este gráfico podría ser:

# Creamos un nuevo factor ordenado de acuerdo a la variable que estamos midiendo
storm$type2 <- reorder(storm$type, storm$wind)
# Creamos el gráfico
ggplot(storm, aes(x = year_f, y = wind, color = type2))  +
  geom_boxplot() + 
  xlab("Año") +
  ylab("Velocidad del viento (mph)")

Un factor, Dos numéricas

En este caso generalizamos el análisis de correlación y el gráfico de dispersión con la inclusión del factor. Consideramos las variables type, wind y pressure. En primer lugar realizamos el estudio descriptivo numérico.

storm %>% 
  group_by(type) %>% 
  summarise(correlación = cor(wind,pressure))

El resultado muestra gran asociación entre la velocidad del viento y la presión atmosférica en todas la categorías salvo para las Depresiones tropicales.

Veamos ahora el gráfico de dispersión conjunto. En primer lugar realizamos un único gráfico marcando con colores los tipos de tormenta

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

Podemos ver como cada punto viene identificado según el tipo de tormenta. Los huracanes en la parte inferior donde se dan las relaciones entre velocidades del viento más altas y presiones atmosféricas más bajas. ¿Qué más podemos decir?

Podemos distinguir cada grupo introduciendo un gráfico matricial

# 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, y = pressure))  +
  geom_point() + 
  xlab("Velocidad del viento (mph)") +
  ylab("Presión atmosférica") +
  facet_grid(. ~ type2)

En este gráfico se aprecia mejor el comportamiento de ambas variables en cada uno de los niveles del factor. Salvo en las depresiones tropicales donde no se aprecia asociación, en el resto de niveles se aprecia un relación de orden inverso. Lo que si podemos ver es que hay observaciones en algunas categorías que podrían corresponder a otras. En la categoría de tormentas tropicales tenemos combinaciones de velocidad y presión que parecen corresponder más a un huracán que a una tormenta tropical. Esto puede ser debido al protocolo de clasificación establecido o a la propia evolución de las tormentas.

Dos factores, Dos numéricas

Este caso es una generalización directa del caso anterior, ya que únicamente debemos añadir una nueva variable categórica. Consideramos las variables type, year_f, wind y pressure. Comenzamos con el análisis numérico:

tabla_cor <- storm %>% 
  group_by(year_f,type) %>% 
  summarise(correlación = cor(wind,pressure))
# Visulaizamos la tabla de resultados de forma óptima
tabla_resumen <- select(tabla_cor,year_f,type,correlación)
spread(tabla_resumen, key = type, value = correlación)

En esta tabla aparecen representados los coeficientes de correlación entre viento y presión para las diferentes combinaciones de niveles de las variables tipo y año. Se aprecian valores muy bajos en todas las combinaciones de la depresión tropical, mientras que en el resto hay asociaciones que pueden resultar interesantes de estudiar posteriormente. En el caso de los huracanes esas asociaciones son muy fuertes ya que muestran valores muy próximos a -1.

En cuanto a los procedimientos gráficos optamos por una combinación de los gráficos que utilizamos en la sección anterior. Representamos el gráfico de dispersión coloreando por tipo de tormenta, y usamos un diagrama matricial por año.

# 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, y = pressure, color = type2))  +
  geom_point() + 
  xlab("Velocidad del viento (mph)") +
  ylab("Presión atmosférica") +
  facet_grid(. ~ year_f)

En todos los años se observa un comportamiento similar del resto de variables, indicando que el año no es una factor que pueda ser considerado como relevante. Si utilizamos la variable mes en su lugar el gráfico resultante es:

# 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, y = pressure, color = type2))  +
  geom_point() + 
  xlab("Velocidad del viento (mph)") +
  ylab("Presión atmosférica") +
  facet_grid(. ~ month_f) +
  theme(axis.text.x = element_text(angle = 90))

En este caso el comportamiento de los meses no es tan parecido. Si bien es cierto que en todas las combinaciones se aprecia una tendencia negativa (sube viento - baja presión), también se puede ver que los huracanes aparecen mayoritariamente en los meses de agosto a octubre. ¿Qué otra información podemos extraer de este gráfico?

Tres factores, Dos numéricas

Es una generalización directa de los dos casos anteriores. Se presenta únicamente el código para obtener los resultados. Añadimos la variable month_f a las del caso anterior.

tabla_cor <- storm %>% 
  group_by(year_f,month_f,type) %>% 
  summarise(correlación = cor(wind,pressure))
# Visulaizamos la tabla de resultados de forma óptima
tabla_resumen <- select(tabla_cor,year_f,month_f,type,correlación)
spread(tabla_resumen, key = type, value = correlación)

¿Qué podemos decir de los resultados obtenidos?

Veamos ahora el gráfico matricial. Seleccionamos los meses de agosto a octubre para poder visualizarlo mejor.

storm_meses <- storm %>%
  filter(month_f == c("August","September","October"))
# Creamos un nuevo factor ordenado de acurdo a la variable que estamos midiendo
storm_meses$type2 <- reorder(storm_meses$type, storm_meses$wind)
# Creamos el gráfico
ggplot(storm_meses, aes(x = wind, y = pressure, color = type2))  +
  geom_point() + 
  xlab("Velocidad del viento (mph)") +
  ylab("Presión atmosférica") +
  facet_grid(year_f ~ month_f) +
  theme(axis.text.x = element_text(angle = 90))

¿Qué podemos decir de este gráfico?

Análisis exploratorio de más de dos variables