1  Descripción de datos

There is nothing better than a picture for making you think of questions you had forgotten to ask.

John W. Tukey

1.1 Introducción

Nuestro objetivo principal en este tema es describir y explorar un conjunto de datos, es decir, presentar la información que contiene de la manera más clara posible, y buscar posibles patrones o relaciones en las variables que lo componen. Necesitaremos examinar la estructura de los datos, detectar los rasgos principales de la distribución de las variables y analizar las relaciones entre ellas usando un conjunto reducido de gráficos y medidas numéricas. Para ello, es imprescindible usar algún software apropiado por lo que un segundo objetivo de este tema es introducir las librerías dplyry ggplot2 del software R, que son especialmente adecuadas para estos propósitos.

En general los conjuntos de datos que vamos a considerar proceden de medir una serie de variables en un subconjunto de individuos (o unidades bajo estudio) de una población. Para que un fichero de datos esté ordenado se deben cumplir los siguientes principios:

• Cada columna corresponde a la información de una variable.
• Cada fila corresponde a la información de un individuo.

Cada fila es una observación, el número de observaciones (filas) es el tamaño muestral (\(n\)) y el número de columnas es la dimensión (\(p\)). Así, un fichero al final toma la forma de una matriz con \(n\) filas y \(p\) columnas. Tradicionalmente la situación más frecuente era disponer de más observaciones que variables, \(n>p\), pero en la actualidad es muy frecuente encontrar conjuntos de datos de alta dimensión para los que \(n<p\). En genética, por ejemplo, en habitual medir los niveles de expresión de miles de genes en unas decenas de individuos.

En términos generales las variables pueden ser:

• variables cualitativas: contienen información sobre cualidades o atributos, que no se expresa numéricamente. Por ejemplo, la variable grupo sanguíneo. En este caso las variables se llaman factores y los valores que toman se llaman niveles.

• variables cuantitativas: contienen información numérica. A su vez, estas variables pueden ser:

  • cuantitativas discretas: si toman un número pequeño de valores, normalmente enteros. Por ejemplo, la variable número de personas que viven en un hogar.
  • cuantitativas continuas: si toman cualquier valor en un intervalo de números reales. Por ejemplo, la variable duración de un componente electrónico hasta que se estropea.

Las variables ordinales mezclan las características de los tipos anteriores. Son variables que toman un número pequeño de valores, numéricos solo de forma imprecisa, entre los cuales se puede definir un orden de manera natural. Por ejemplo, cuando se aplica un tratamiento en el que se administra un medicamento en forma de dosis baja, dosis media y dosis alta.

Para describir un conjunto de datos es aconsejable usar en primer lugar técnicas gráficas y posteriormente resúmenes numéricos. Es importante tener en cuenta el tipo de cada variable a la hora de elegir la manera de describirla.

1.2 Distribución de una variable

La distribución de una variable en un conjunto de datos viene determinada por los valores que toma esa variable y la frecuencia con la que los toma. Un primer objetivo de explorar un conjunto de datos es describir los rasgos principales de las distribuciones de las variables que contiene. Si consideramos que los datos \(x_1,\ldots,x_n\) corresponden a \(n\) realizaciones independientes de una variable aleatoria entonces la distribución de los datos es una aproximación de la distribución de la variable aleatoria. Diremos mucho más sobre esto en los capítulos siguientes puesto que, en sentido amplio, un objetivo genérico de la estadística es inferir las propiedades de la distribución de la variable aleatoria a partir de la información parcial que proporciona la distribución de los datos.

La frecuencia absoluta de un valor (o de un intervalo) es el número de observaciones para las que la variable toma ese valor (o pertenece a ese intervalo). Las frecuencias absolutas se denotarán \(N_1,\ldots,N_k\), donde \(k\) es el número de valores diferentes que aparecen en los datos.

La frecuencia relativa es igual a la frecuencia absoluta dividida por el número de datos \(n\). Siempre toma valores entre cero y uno. Las denotaremos \(f_1,\ldots,f_k\).

En ocasiones nos encontraremos con datos agrupados en \(k\) intervalos o clases \(A_1,\dots,A_k\). Los valores \(m_1,\dots,m_k\) que representan cada clase (generalmente los centros de los intervalos) se suelen llamar marcas de clase.

De la distribución de una variable cuantitativa nos interesa conocer:

• Su posición: en torno a qué valor central toma valores la variable.
• Su dispersión: el grado de concentración de los valores que toma la variable alrededor de su posición central.
• Su forma: por ejemplo, la simetría, es decir, si los valores se reparten de la misma forma a uno y otro lado del centro.

Cuestiones

1. ¿Qué diferencia hay entre \(n\) y \(k\)? ¿Para qué tipo de variables podemos esperar que \(n\approx k\)? ¿Cuándo serán muy diferentes?
2. ¿Cuánto vale \(N_1+\cdots + N_k\)?
3. ¿Cuánto vale \(f_1+\cdots + f_k\)?
4. Piensa en dos conjuntos de 5 datos que tengan

• la misma posición y distinta dispersión,
• la misma dispersión y distinta posición.

1.3 La estructura fundamental de datos en R

Antes de comenzar a usar los gráficos y las medidas numéricas más importantes para describir una distribución, es necesario aprender a manipular un conjunto de datos con R. La estructura básica para organizar los datos es un data frame, que corresponde a un fichero de datos con \(n\) observaciones y \(p\) variables. Cada una de las variables es un vector (que puede contener información numérica o cualitativa), pero todas las variables deben tener la misma longitud. Por lo tanto un data frame es una matriz de datos, pero se diferencia de las matrices en que las columnas pueden ser de distinta naturaleza (números, palabras, fechas, etc.) En este sentido, también es posible entender un data frame como una lista de variables, una lista restringida a que todos sus elementos (las variables) tienen la misma longitud (el número de observaciones).

Para crear un data frame con R se usa el comando data.frame. En este ejemplo creamos un fichero llamado df con dos variables, una cuantitativa y otra cualitativa:

x <- c(1, 2, 3, 4)  # variable cuantitativa
y <- c("a", "b", "c", "d")  # variable cualitativa

datos <- data.frame(var1 = x, var2 = y)  # se crea el data frame
glimpse(datos)   # para ver la estructura del fichero
#> Rows: 4
#> Columns: 2
#> $ var1 <dbl> 1, 2, 3, 4
#> $ var2 <chr> "a", "b", "c", "d"

Habitualmente no tendremos que introducir nosotros los datos, sino que estarán ya disponibles en diferentes formatos (texto, excel, etc.) y tendremos que crear el data frame a partir de cada formato. Por ejemplo, el comando read.table es el básico para leer datos de un fichero de texto, aunque hay multitud de comandos y paquetes para este propósito.

Ejemplo

El agua de los ríos contiene pequeñas concentraciones de mercurio que se pueden ir acumulando en los tejidos de los peces. Se ha realizado un estudio en los ríos Wacamaw y Lumber en Carolina del Norte (EE.UU.), analizando la cantidad de mercurio que contenían 171 ejemplares capturados de una cierta especie de peces. El siguiente código lee el fichero de texto (situado en una página web) que contiene los datos y lo convierte en un data.frame. La opción header = TRUE sirve para informar de que la primera fila del fichero de texto ya contiene los nombres de las variables (por defecto, se asume que no es así):

enlace <- "http://matematicas.uam.es/~joser.berrendero/datos/mercurio.txt"
mercurio <- read.table(enlace, header = TRUE) 

glimpse(mercurio)
#> Rows: 171
#> Columns: 5
#> $ RIO      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
#> $ ESTACION <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
#> $ LONG     <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4, 4~
#> $ PESO     <dbl> 1616, 1862, 2855, 1199, 1320, 1225, 870, 1455, 1220, 1033, 33~
#> $ CONC     <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50, 0~

El fichero consta de cinco variables: el río en el que se ha capturado cada pez (0 corresponde al río Lumber y 1 corresponde al río Wacamaw), la estación o punto de cada río en el que se ha capturado el pez (codificada de 0 a 15), la longitud (cm), el peso (g) y la concentración de mercurio en los tejidos de cada pez (ppm).

En este ejemplo, ¿de qué tipo es cada una de las variables? ¿Cuántas observaciones hay? ¿De qué dimensión?

1.4 Manejo básico de los datos

Vamos a usar el paquete dplyr para manejar los datos. Está estructurado en una serie de verbos que corresponden a las manipulaciones más habituales. Incorpora una sintaxis en cadena (pipes) que sigue este esquema:

nuevo_data_frame <- antiguo_data_frame %>% 
  verbo1("argumentos") %>% 
  verbo2("argumentos") %>% 
  ...

El operador %>% sirve para encadenar secuencialmente las acciones que se van llevando a cabo sobre los datos. Podemos traducirlo por la palabra entonces. Lo que aparece a la izquierda de este operador es el primer argumento del comando que aparece a su derecha. En el ejemplo siguiente se compara la sintaxis en cadena con la habitual sintaxis anidada:

# Calcula la distancia euclídea entre dos vectores x1 y x2
x1 <- 1:6
x2 <- 7:12

# En cadena
(x1 - x2)^2 %>%
  sum() %>%
  sqrt() 
#> [1] 14.69694

# Anidada
sqrt(sum((x1-x2)^2))
#> [1] 14.69694

Normalmente la sintaxis en cadena proporciona un código mucho más fácil de leer y de entender. El operador %>%está implementado en el paquete magrittr y se puede usar una vez cargado tidyverse. Las versiones más recientes de R incorporan un operador nativo |>, que se usa de la misma forma:

(x1 - x2)^2 |>
  sum() |>
  sqrt()
#> [1] 14.69694

Los principales verbos son: mutate, count, filter, select, arrange y summarise. Veamos para qué sirve cada uno de ellos.

1.4.1 Mutate

Para modificar las variables existentes o crear otras nuevas que resultan de transformarlas.

Lo primero que vamos a hacer es recodificar la variable RIO, de forma que aparezcan los nombres de los ríos en lugar de los códigos 0 y 1, y vamos a definir tanto esta variable como la variable ESTACIÓN como variables cualitativas (factores). El resultado lo enviamos a un data frame que se llama igual que el original y por lo tanto lo sustituye:

mercurio <- mercurio %>% 
  mutate(RIO = recode(RIO, "0" = "Lumber", "1" = "Wacamaw")) %>% 
  mutate(RIO = factor(RIO)) %>% 
  mutate(ESTACION = factor(ESTACION))

glimpse(mercurio)
#> Rows: 171
#> Columns: 5
#> $ RIO      <fct> Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumbe~
#> $ ESTACION <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
#> $ LONG     <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4, 4~
#> $ PESO     <dbl> 1616, 1862, 2855, 1199, 1320, 1225, 870, 1455, 1220, 1033, 33~
#> $ CONC     <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50, 0~

1.4.2 Count

Para contar, es decir, obtener las frecuencias absolutas de cada valor o combinación de valores.

mercurio %>% 
  count(RIO)
#>       RIO  n
#> 1  Lumber 73
#> 2 Wacamaw 98

El comando anterior nos da el número de peces capturados en cada río. Es importante observar que la salida tiene también la estructura de data frame, con una variable (llamada RIO) con los dos niveles y otra (llamada n) con las dos frecuencias. Si usamos el comando con dos variables en lugar de una nos da las frecuencias de cada combinación de niveles de las dos variables:

mercurio %>% 
  count(RIO, ESTACION)
#>        RIO ESTACION  n
#> 1   Lumber        0 10
#> 2   Lumber        1 19
#> 3   Lumber        2 13
#> 4   Lumber        3 10
#> 5   Lumber        4 14
#> 6   Lumber        5  3
#> 7   Lumber        6  4
#> 8  Wacamaw        7 10
#> 9  Wacamaw        8 14
#> 10 Wacamaw        9  3
#> 11 Wacamaw       10 21
#> 12 Wacamaw       11  9
#> 13 Wacamaw       12 13
#> 14 Wacamaw       13  4
#> 15 Wacamaw       14 13
#> 16 Wacamaw       15 11

Ahora el data frame resultante tiene tres variables, dos para los niveles de las variables y otra para la frecuencia de cada combinación de niveles.

1.4.3 Filter

Para filtrar las observaciones y quedarnos solo con las que cumplen cierta condición.

Por ejemplo, nos quedamos solo con las observaciones que corresponden al río Lumber:

mercurio_lumber <- mercurio %>% 
  filter(RIO == "Lumber")

glimpse(mercurio_lumber)
#> Rows: 73
#> Columns: 5
#> $ RIO      <fct> Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumbe~
#> $ ESTACION <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
#> $ LONG     <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4, 4~
#> $ PESO     <dbl> 1616, 1862, 2855, 1199, 1320, 1225, 870, 1455, 1220, 1033, 33~
#> $ CONC     <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50, 0~

En este otro ejemplo retenemos solo las observaciones del río Lumber tales que la longitud es mayor que 45:

datos <- mercurio %>% 
  filter(RIO == "Lumber" & LONG > 45)

glimpse(mercurio_lumber)
#> Rows: 73
#> Columns: 5
#> $ RIO      <fct> Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumbe~
#> $ ESTACION <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
#> $ LONG     <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4, 4~
#> $ PESO     <dbl> 1616, 1862, 2855, 1199, 1320, 1225, 870, 1455, 1220, 1033, 33~
#> $ CONC     <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50, 0~

1.4.4 Select

Para seleccionar un subconjunto de las variables.

Si queremos retener solo las variables correspondientes a la longitud y la concentración de mercurio:

datos <- mercurio %>% 
  select(LONG, CONC) 

glimpse(datos)
#> Rows: 171
#> Columns: 2
#> $ LONG <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4, 47.7,~
#> $ CONC <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50, 0.80,~

Existe una sintaxis que se puede consultar en este enlace para seleccionar variables de manera eficiente. Por ejemplo, es posible seleccionar todas las variables numéricas, todas las que son factores, todas aquellas cuyo nombre incluya una cierta palabra, etc.

1.4.5 Arrange

Para ordenar los datos.

En el siguiente ejemplo seleccionamos los peces del río Lumber cuya longitud está por encima de la media de las longitudes del río y los ordenamos de menor a mayor concentración de mercurio:

datos <- mercurio_lumber %>% 
  filter(LONG > mean(LONG)) %>% 
  arrange(CONC)

glimpse(datos)
#> Rows: 35
#> Columns: 5
#> $ RIO      <fct> Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumbe~
#> $ ESTACION <fct> 1, 0, 0, 1, 0, 4, 1, 1, 0, 4, 1, 1, 0, 3, 3, 3, 3, 3, 5, 0, 0~
#> $ LONG     <dbl> 45.1, 40.4, 43.8, 43.5, 44.7, 41.0, 47.4, 44.2, 45.2, 44.0, 4~
#> $ PESO     <dbl> 2920, 1033, 1225, 2674, 1320, 1042, 3675, 1378, 1199, 1455, 3~
#> $ CONC     <dbl> 0.34, 0.50, 0.51, 0.54, 0.56, 0.66, 0.69, 0.72, 0.73, 0.73, 0~

1.4.6 Summarise

Para calcular medidas descriptivas básicas.

Veremos ejemplos de uso de este verbo más adelante, en la sección dedicada a las medidas numéricas.

1.5 Descripción gráfica de una variable

1.5.1 Nociones básicas de ggplot2

Para visualizar un conjunto de datos vamos a usar el paquete ggplot2 de R. El esquema esencial para representar un gráfico con ggplot2 es el siguiente:

ggplot(data = <DATA>) + 
  <GEOM_FUNCTION>(mapping = aes(<MAPPINGS>))

El comando ggplot(data = <DATA>) crea un sistema de coordenadas vacío. El argumento data determina el data frame donde están los datos que se quieren representar. Es importante saber que los datos siempre tienen que tener estructura de data frame. Esta rigidez aporta sin embargo claridad sobre lo que se está haciendo y facilita la reproducibilidad de los resultados.

Al sistema de coordenadas vacío creado se le van añadiendo capas usando el signo +.

El comando <GEOM_FUNCTION>() añade una capa al gráfico en la que se determina el tipo de elementos que se van a representar: puntos, líneas, histogramas, diagramas de cajas, etc. Hay un catálogo muy amplio de geoms que se puede consultar en la documentación del paquete.

El argumento mapping = aes(<MAPPINGS>) determina los llamados aesthetics, una lista en la que se asignan las variables del fichero a distintos aspectos del gráfico. Por ejemplo, casi siempre hay que asignar el eje de abscisas a una variable y el eje de ordenadas a otra. Cada geom admite una lista diferente de aesthetics que se puede consultar en la documentación correspondiente.

Se pueden añadir todas las capas que resulten necesarias, correspondientes a varios geoms o a otros aspectos del gráfico (como por ejemplo título y etiquetas, o alguna transformación de las variables originales). Iremos viendo ejemplos de uso para ir comprendiendo como funcionan estas capas.

Las posibilidades de visualización deggplot2 se potencian a través de extensiones programadas por los usuarios. Un catálogo de las extensiones disponibles se puede encontrar en esta página.

1.5.2 Diagrama de barras

Para representar gráficamente la distribución de una variable cualitativa (o a veces tambien de una variable cuantitativa discreta) lo más sencillo es usar un gráfico de barras. Sobre cada uno de los valores o niveles de la variable se representa una barra cuya longitud es proporcional a la frecuencia de veces con que se repite ese valor. A continuación vemos el diagrama de barras para la tabla de frecuencias correspondiente a los peces capturados en cada estación:

# tabla de frecuencias
mercurio %>% 
  count(RIO)
#>       RIO  n
#> 1  Lumber 73
#> 2 Wacamaw 98

# gráfico de barras
ggplot(mercurio) +
  geom_bar(aes(x = ESTACION, fill = RIO),
           width = 0.5) +
  labs(x = "Estación", y = "Número de capturas")

Hemos hecho depender el color de las barras del río en el que se ha hecho la captura (asignando RIO a fill en la lista de aesthetics). Vemos que las estaciones 0 a 6 están situadas en el río Lumber mientras que las estaciones 7 a 15 están en el Wacamaw. Obsérvese también el uso en la última capa del comando labs() para cambiar las etiquetas de los ejes.

Con el mismo propósito que los diagramas de barras a veces se usan los gráficos de tartas pero, a pesar de su uso generalizado, no se recomiendan porque son más difíciles de interpretar que los de barras. Se atribuye a John Tukey la frase:

There is no data that can be displayed in a pie chart, that cannot be displayed BETTER in some other type of chart.

¿Qué ocurre si calculamos la tabla de frecuencias o representamos el diagrama de barras de una variable continua?

1.5.3 Histogramas

Para construir un histograma:

1. Se divide el rango de los datos en un número adecuado de intervalos.
2. Sobre cada intervalo se dibuja un rectángulo cuya área es proporcional a la frecuencia (relativa o absoluta) de datos en el intervalo.

Este es el histograma de la variable correspondiente a la concentración de mercurio:

ggplot(mercurio) +
  geom_histogram(aes(x = CONC, y = after_stat(density)),  # after_stat(density) normaliza para área = 1
                 bins = 10,                       # número de intervalos
                 col = 'black',                   # color del borde
                 fill = 'olivedrab4') +           # color del interior
  labs(x = "Concentración (ppm)", y = NULL)        # cambia etiquetas

La asignación y = after_stat(density) es aconsejable para fijar la escala del gráfico de forma que la suma de todas las áreas de los rectángulos del histograma sea igual a 1 y por lo tanto el histograma sea comparable con una función de densidad. El otro argumento relevante es bins, que permite determinar el número de intervalos.

Cuestiones

Un histograma no es lo mismo que un diagrama de barras. Enumera sus semejanzas y diferencias.

Algunos aspectos a tener en cuenta cuando miramos un histograma:

• Una moda es un máximo relativo del histograma. Corresponde a un valor de alta frecuencia rodeado de valores de frecuencias más bajas. Nos tenemos que fijar en el número de modas. La presencia de varias modas suele indicar que hay varias poblaciones mezcladas (por ejemplo a causa de una variable latente que no estamos teniendo en cuenta y que influye en la variable del histograma). También hay que distinguir en la medida de lo posible las modas que corresponden a la población de la que proceden los datos de aquellos máximos que corresponden a ruido en la muestra disponible.
• Si hay datos atípicos o no en la muestra. Por un lado, la presencia de atípicos puede distorsionar otros cálculos que hagamos posteriormente. Por otro lado, a veces un dato atípico puede llevarnos a descubrir algún fenómeno en la población del que no éramos conscientes.
• Si la distribución es o no simétrica y qué tipo de simetría tiene. Esto tiene interés para interpretar correctamente las medidas numéricas y más adelante para ver si los datos se ajustan a algunos modelos teóricos. En el ejemplo del histograma anterior la distribución presenta asimetría a la derecha. Las variables económicas (ingresos) pueden ser muy asimétricas a la derecha. Menos frecuente es la asimetría a la izquierda aunque también hay ejemplos (la esperanza de vida medida en diferentes países). La distribución de medidas de una misma magnitud (por ejemplo, la temperatura en un lugar medida con varios termómetros distintos) suele tener la típica forma simétrica de campana, que sugiere una distribución normal como modelo teórico apropiado.
• Cuál es la magnitud en torno a la cual se distribuyen los datos. Si hay mucha dispersión en torno a este centro.

La forma de una histograma depende del número de clases, tal y como vemos en la figura siguiente:

Tenemos que elegir el número de intervalos que mejor creemos que representaría los rasgos principales de la distribución en la población de la que proceden los datos. Una elección de un número demasiado bajo de intervalos puede ocultar rasgos importantes. Una elección de un número demasiado alto nos lleva a un histograma que refleja los datos concretos que tenemos pero no la población de la que proceden (se dice que el procedimiento sobreajusta los datos).

1.5.4 Estimadores de núcleo

Podemos considerar los histogramas como aproximaciones o estimaciones de la función de densidad de la que proceden los datos. Sin embargo, un histograma no es del todo satisfactorio ya que su aspecto (son funciones constantes a trozos) no corresponde al de las funciones de densidad más habituales. En este apartado consideramos una clase de estimadores un poco más sofisticados.

Sabemos que las áreas bajo la función de densidad corresponden a las probabilidades de que la variable tome valores en el correspondiente intervalo. Por lo tanto, si \(h\approx 0\) y \(f\) es la función de densidad de la variable aleatoria \(X\) se tiene la siguiente aproximación: \[\mbox{P}(x-h\leq X\leq x+h) = \int_{x-h}^{x+h} f(t) dt \approx 2h f(x).\] Usando esta observación construimos un estimador de la función de densidad prefijando \(h>0\) pequeño y contando la proporción de datos que caen dentro del intervalo \((x-h, x+h)\). El estimador resultante es \[\hat{f}(x) = \frac{1}{2h} \frac{\#\{i:\, |x-x_i| \leq h\}}{n} = \frac{1}{2hn} \#\left\{i:\, \frac{|x-x_i|}{h} \leq 1\right\}.\] Si definimos la función \(K(x) = (1/2) \mathbb{I}_{\{|x|\leq 1\}}\), donde \(\mathbb{I}_A\) denota la función indicatriz sobre \(A\), que vale 1 en \(A\) y 0 fuera de \(A\), entonces podemos reescribir el estimador como: \[\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^nK\left(\frac{x-x_i}{h}\right)= \frac{1}{n} \sum_{i=1}^n \frac{1}{h}K\left(\frac{x-x_i}{h}\right).\] Fijado un valor de \(x\), comprobamos si pertenece al intervalo centrado en \(x_i\) y de radio \(h\). Asignamos un valor \(1/(2h)\) cada vez que esto ocurre y cero en caso contrario, y entonces promediamos los valores. Nótese que la función \(K\) que hemos definido corresponde a la función de densidad de una variable aleatoria (v.a.) uniforme en el intervalo \((-1,1)\).

Si sustituimos la función indicatriz por una función más suave obtendremos estimadores que se asemejan más a una función de densidad típica. Así se obtienen los estimadores de núcleo de la densidad: \[\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^n K\left(\frac{x-X_i}{h}\right),\] donde \(h>0\) es el parámetro de suavizado y \(K\) es el núcleo (que suele verificar \(K\geq 0\) y \(\int K = 1\), es decir, que es una función de densidad). Si, por ejemplo, \(K(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}\) es un núcleo gaussiano (la función de densidad de una v.a. normal estándar) entonces los valores \[ \frac{1}{h}K\left(\frac{x-x_i}{h}\right) \] corresponderán a las funciones de densidad de v.a. con distribución normal de media \(x_i\) y desviación típica \(h\). El estimador núcleo es entonces el promedio de los valores de estas densidades en cada punto \(x\). En la siguiente figura se han representado las densidades \(h^{-1}K((x-x_i)/h)\) divididas por \(n\) de manera que la suma es el estimador núcleo (en azul):

Si el núcleo \(K\) es una función de densidad, también lo es \(\hat{f}\) (ya que es muy fácil comprobar que \(\hat{f}\geq 0\) y \(\int \hat{f}(x) dx = 1\)).

Usamos el comando geom_density para representar los estimadores de núcleo. El comando geom_density usa por defecto un núcleo gaussiano:

ggplot(mercurio) +
  geom_density(aes(x = CONC),
                 fill = 'olivedrab4') +
  labs(x = "Concentración (ppm)", y = " ")

Si usamos un núcleo rectangular el resultado es menos suave, más parecido al histograma, pero los rasgos básicos del estimador se mantienen:

ggplot(mercurio) +
  geom_density(aes(x = CONC),
               fill = 'olivedrab4',
               kernel = 'rectangular') +
  labs(x = "Concentración (ppm)", y = " ")

La selección del parámetro de suavizado \(h\) es un problema difícil e importante ya que el valor de este parámetro sí determina drásticamente las propiedades del estimador. Su papel es análogo al del número de intervalos en un histograma. En términos generales puede comprobarse que si queremos buenas propiedades del estimador (lo que significa que -en un cierto sentido que se puede precisar- esté cerca de la verdadera función de densidad) debe cumplirse que \(h=h_n\) sea del orden de \(n^{-1/5}\), esto quiere decir que \(h_n\) debe decrecer con el tamaño muestral, pero bastante lentamente. Un poco más adelante comentaremos la opción por defecto que usa geom_density para elegir el parámetro de suavizado. La siguiente figura muestra cómo cambia el estimador para diferentes valores de \(h\):

Como técnica descriptiva, una ventaja sencilla pero relevante de los estimadores de núcleo respecto a los histogramas es que son más apropiados para superponerlos en el mismo gráfico, de manera que se puedan comparar más fácilmente dos distribuciones:

ggplot(mercurio) +
  geom_density(aes(x = CONC, linetype = RIO)) +
  labs(x = "Concentración (ppm)", y = " ")

Los gráficos de cordilleras son una forma más sofisticada de comparar estimadores de núcleo. Están implementados en el paquete ggridges:

# Requiere el paquete ggridges

ggplot(mercurio) +
  geom_density_ridges(aes(x = CONC, y = RIO, fill = RIO)) +
  labs(x = "Concentración (ppm)", y = "Río")

1.6 Descripción numérica de una variable

A la hora de describir numéricamente diversos aspectos de una distribución de probabilidad, debemos distinguir entre medidas de posición, dispersión y forma (por ejemplo, asimetría o simetría).

1.6.1 Posición: media, mediana y percentiles

La media de un conjunto de números \(x_1,\ldots,x_n\) se define como \(\bar{x}=(x_1+\cdots + x_n)/n\). En general esta medida es un centro en el sentido de que las desviaciones positivas se compensan con las negativas y la suma de las desviaciones de los datos a la media es igual a cero: \[(x_1-\bar{x}) + \cdots + (x_n -\bar{x}) = 0.\]

x <- c(1, 2, 5, 8)
mean(x)
#> [1] 4
x - mean(x)
#> [1] -3 -2  1  4
sum(x - mean(x))
#> [1] 0

Como consecuencia de esta propiedad de las desviaciones, la media se puede interpretar como el centro de gravedad de un histograma, en el sentido de que sería el lugar en el que apoyar el histograma para que se mantuviera en equilibrio si los pesos de los rectángulos que lo forman fueran proporcionales a sus áreas.

Otra consecuencia es que \(\bar{x}\) es el valor que minimiza la suma de desviaciones al cuadrado, es decir, si \[f(t) := (x_1-t)^2 + \cdots + (x_n-t)^2, \tag{1.1}\] entonces \(f'(\bar{x})=0\), donde \(f'\) es la derivada, y por tanto \(\bar{x}\) es el valor donde \(f\) alcanza el mínimo.

Hay dos situaciones en las que la media puede dar una idea equivocada de la posición:

• Si la distribución es muy asimétrica (aquí no está muy claro cómo interpretar la media)
• Si hay datos atípicos (distorsionan su valor)

Una medida de posición que se comporta mejor respecto a estos dos problemas es la mediana. Para calcularla:

• Se ordenan los datos de menor a mayor.
• Si el número de datos es impar, la mediana es el dato que ocupa la posición central. Si es par, la mediana es la media de los dos datos centrales (se puede entender también que todos los valores entre los dos datos centrales son medianas, con lo que estaríamos eligiendo un representante razonable de ellos).

Por ejemplo las medianas de 1, 2, 5, 8 serían los valores del intervalo (2, 5) o directamente 3.5 si queremos dar un valor único. La mediana de 1, 2, 7, 4, 0 es 2.

Como ilustra el meme de la figura 1.1, la mediana es más robusta que la media (menos sensible a datos atípicos). Por otra parte, hace un uso menos eficiente de la información contenida en los datos (más sobre esto en el tema siguiente).

Figura 1.1: Los datos atípicos influyen mucho más en la media que en la mediana (fuente: cuenta de Twitter de annaegalite)

Puede demostrarse que la mediana es el valor de \(t\) en el que se minimiza la función \[g(t):= |x_1-t| +\cdots + |x_n-t|.\] Compara esta función con la función definida en (1.1), cuyo mínimo se alcanza en la media.

Si una distribución es simétrica, la media y la mediana coinciden aproximadamente. Si es asimétrica a la derecha, como en el caso de muchas variables económicas -véase la figura 1.2 tomada de la página del INE)- la media suele ser mayor que la mediana. Si es asimétrica a la izquierda la media suele ser menor que la mediana.

Figura 1.2: Media, mediana y moda de los salarios en 2018 (fuente: Instituto Nacional de Estadística)

El cuantil \(p\) (o percentil \(100p\)) es el valor que deja a su izquierda una proporción \(p\) (o porcentaje \(100p\)) de los datos.
Al percentil 25 se le llama también primer cuartil o cuartil inferior (notación: \(Q_1\)). El segundo cuartil es la mediana (que también es el percentil 50). Al percentil 75 se le suele llamar tercer cuartil o cuartil superior (notación: \(Q_3\)). En la figura 1.3 se muestran algunos percentiles correspondientes a salarios en España en 2017, tomados de la encuesta de estructura salarial del INE:

Figura 1.3: Datos de salarios de 2017 (fuente: Instituto Nacional de Estadística)

Cuestiones

Calcula el valor de \(x\) y de \(y\) en cada frase:

• La media y la mediana de los salarios en España en 2017 fueron \(x\) e \(y\), respectivamente.
• El rango intercuartílico de los salarios en España en 2017 es \(x\).
• Un 10% de las mujeres ganaba más de \(x\) euros.
• Un \(x\%\) de lo hombres ganaba entre \(11393.30\) y \(45813.46\) euros.
• Un 90% de las mujeres ganaba más de \(x\) euros.

1.6.2 El verbo summarise de dplyr

Para calcular las medidas descriptivas más básicas se puede usar el verbo summarise de dplyr (usualmente en combinación con group_by). El comando group_by divide los datos según los niveles de un factor de manera que todos los cálculos posteriores se realizan por separado para cada nivel. En el siguiente ejemplo se calculan la media, \(Q_1\), la mediana y \(Q_3\) para la variable correspondiente a la concentración de mercurio y para cada río por separado:

mercurio %>% 
  group_by(RIO) %>% 
  summarise(media_hg = mean(CONC),
            Q1_hg = quantile(CONC, probs = 0.25),
            mediana_hg = median(CONC),
            Q3_hg = quantile(CONC, probs = 0.75))
#> # A tibble: 2 x 5
#>   RIO     media_hg Q1_hg mediana_hg Q3_hg
#>   <fct>      <dbl> <dbl>      <dbl> <dbl>
#> 1 Lumber      1.08 0.57        0.91   1.5
#> 2 Wacamaw     1.28 0.618       1.05   1.8

Comparando las filas del resultado vemos que las cantidades de mercurio en los peces del Wacamaw tienden a ser mayores que en el Lumber (cosa que ya habíamos podido apreciar en los gráficos de los estimadores de las densidades).

El comando across() sirve para trabajar con muchas variables y medidas descriptivas a la vez. Es un comando muy potente que toma dos argumentos: una lista de variables y una lista de medidas numéricas. A continuación vemos un ejemplo en el que se calculan, para todas las variables numéricas del fichero y para cada río, la media y la mediana:

mercurio %>% 
  group_by(RIO) %>% 
  summarise(across(LONG:CONC, list(media = mean, mediana = median)))
#> # A tibble: 2 x 7
#>   RIO     LONG_media LONG_mediana PESO_media PESO_mediana CONC_media CONC_medi~1
#>   <fct>        <dbl>        <dbl>      <dbl>        <dbl>      <dbl>       <dbl>
#> 1 Lumber        39.4         39        1197.          898       1.08        0.91
#> 2 Wacamaw       40.4         39.2      1111.          860       1.28        1.05
#> # ... with abbreviated variable name 1: CONC_mediana

En este enlace se puede encontrar más información sobre las posibilidades de across() para operar sobre un subconjunto de las variables del fichero simultáneamente.

1.6.3 Dispersión: varianza, desviación típica y rango intercuartílico

Estas son (sobre todo la segunda) las medidas de dispersión más utilizadas.

La varianza es básicamente el promedio de las desviaciones al cuadrado de los datos a su media. Al calcular la varianza estamos midiendo si las observaciones tienden a estar cerca o lejos de la media.

x <- c(1, 2, 3, 5, 7, 0)
x - mean(x)   # desviaciones
#> [1] -2 -1  0  2  4 -3
(x - mean(x))^2   # desviaciones al cuadrado
#> [1]  4  1  0  4 16  9
mean((x-mean(x))^2) # promedio de desviaciones al cuadrado
#> [1] 5.666667

Como sabemos que la suma de desviaciones siempre es cero, a partir de todas ellas menos una podemos deducir la que falta. Es decir, solo disponemos de \(n-1\) desviaciones independientes. Como consecuencia es más correcto dividir por \(n-1\) que por \(n\). Por ello, en la mayoría de los libros la varianza se define de la siguiente forma (algunos autores llaman a esta cantidad cuasivarianza):

\[s^2 = \frac{(x_1-\bar{x})^2+\cdots +(x_n-\bar{x})^2}{n-1}.\]

Cuestiones

Dadas dos variables aleatorias independientes e idénticamente distribuidas \(X_1\) y \(X_2\) con media común \(\mu\) y varianza común \(\sigma^2\),

• comprueba que \(\mbox{E}[(X_1-\mu)^2 + (X_2-\mu)^2]=2\sigma^2\),
• comprueba que \(\mbox{E}[(X_1-\bar{X})^2 + (X_2-\bar{X})^2]=\sigma^2\).

La cuestión anterior apunta a que, si nuestro objetivo es aproximar \(\sigma^2\) y disponemos de dos datos, debemos dividir por dos las desviaciones \((X_i-\mu)^2\); pero si en lugar de usar \(\mu\) usamos \(\bar{X}\) (como tendremos que hacer siempre en la práctica) no hay que dividir por dos sino por uno (es decir, \(n-1\) en este caso).

Una fórmula alternativa (ejercicio) para \(S^2\) es: \[s^2 = \frac{n}{n-1}\left(\frac{x_1^2+\cdots + x_n^2}{n} - \bar{x}^2 \right).\]

La desviación típica es la raíz cuadrada de \(s^2\). Al tomar la raíz se compensa el hecho de que habíamos considerado las desviaciones al cuadrado de forma que se mide la dispersión en las mismas unidades que los datos: \[s =\sqrt{\frac{(x_1-\bar{x})^2+\cdots +(x_n-\bar{x})^2}{n-1}}.\]

Una medida adimensional relacionada es el coeficiente de variación: \[\mbox{CV} = \frac{s}{|\bar{x}|}\]

Una medida de dispersión menos sensible a datos atípicos que \(s\) es el rango intercuartílico, que se define como \(\mbox{RI}= Q_3-Q_1\).

Para los datos de los peces, las medidas de dispersión que hemos definido se pueden calcular de la forma siguiente:

mercurio %>% 
  group_by(RIO) %>% 
  summarise(var_hg = var(CONC),
            sd_hg = sd(CONC),
            RI_hg = quantile(CONC, probs = 0.75) - quantile(CONC, probs = 0.25))
#> # A tibble: 2 x 4
#>   RIO     var_hg sd_hg RI_hg
#>   <fct>    <dbl> <dbl> <dbl>
#> 1 Lumber   0.421 0.649  0.93
#> 2 Wacamaw  0.687 0.829  1.18

El rango intercuartílico interviene en el criterio automático para seleccionar el parámetro de suavizado de los estimadores de núcleo. De hecho, la regla que usa por defecto geom_density es \[h_n = 0.9 \hat{\sigma} n^{-1/5}\] donde \(\hat{\sigma}=\min\{s,\mbox{RI}/1.34\}\) es una estimación de la desviación típica de la población con un cierto grado de robustez.

Cuestiones

1. ¿Por qué crees que se divide el rango intercuartílico por 1.34 en la regla que usa por defecto geom_density para seleccionar el parámetro de suavizado? (Pista: calcula el rango intercuartílico aproximado de datos cuya distribución sea \(\mbox{N}(\mu,\sigma^2)\)).
2. Siempre se cumple \(s^2\geq 0\). Da un ejemplo de un conjunto de datos tal que \(s^2=0\).
3. Dado un conjunto de observaciones medidas en kg, supongamos que cambiamos las unidades y las pasamos a gramos (es decir, multiplicamos por mil). Determina si son verdaderas o falsas las siguientes afirmaciones:
   • Tanto la media como la mediana de los nuevos datos se multiplican también por mil.
   • La varianza se multiplica también por mil.
4. ¿Cómo cambiaría la desviación típica en el caso anterior?
5. Ahora sumamos 100 a todos los datos. Determina si son verdaderas o falsas las siguientes afirmaciones:
   • Los cuartiles no cambian.
   • El rango intercuartílico no cambia.
   • La desviación típica no cambia.

1.6.4 Estandarización

Dado un conjunto de datos \(x_1,\ldots, x_n\) con media \(\bar{x}\) y desviación típica \(s\), estandarizarlos consiste en restar a cada uno de ellos la media y dividir el resultado entre \(s\). Por tanto el dato \(x_i\) estandarizado es \[z_i = \frac{x_i - \bar{x}}{s}.\]

Cuestiones

1. ¿Cómo cambian los datos estandarizados si se cambia la unidad de medida de la variable a la que corresponden?
2. Dado un conjunto de datos estandarizados arbitrario \(z_1,\ldots,z_n\), calcula su media \(\bar{z}\) y su desviación típica \(s_z\).

El valor \(|z_i|\) puede interpretarse como la distancia entre \(x_i\) y \(\bar{x}\) medida en desviaciones típicas. Si para un dato tenemos \(z_i=2\) podemos decir que \(x_i\) dista de \(\bar{x}\) dos desviaciones típicas. De hecho la desviación típica es una medida adecuada de distancia en estadística puesto que un valor \(|x_i-\bar{x}|\) resulta más significativo cuanto más concentrados están los datos, es decir, cuanto menor es \(s\). Si todos los datos son muy parecidos, una pequeña distancia a la media es estar muy lejos. Por ello, al estandarizar podemos detectar posibles valores atípicos, aquellos para los que \(|z_i|>c\) para un valor adecuado \(c\). Si, por ejemplo, suponemos que los valores proceden de una distribución normal, ¿qué valor de \(c\) te parece adecuado?

Estandarizar también es útil para comparar los valores de una variable en diferentes poblaciones. Por ejemplo, podemos comparar notas de estudiantes de carreras cuya nota media es muy diferente.

1.6.5 Diagramas de cajas

Las partes de un diagrama de cajas se describen en la figura 1.4 tomada de la web towardsdatascience

Figura 1.4: Definición de un diagrama de cajas. (Fuente: página web towardsdatascience.)

Los diagramas de cajas permiten también detectar atípicos. Observa que un dato se marca o no como atípico dependiendo de una medida de dispersión, el rango intercuartílico. Si este es pequeño, una ligera desviación de la mediana puede ser ya sufiente como para marcar un dato como posible atípico.

La principal utilidad de un diagrama de cajas es comparar las características de varios conjuntos de datos. Para representarlo con R se usa geom_boxplot. Por ejemplo, para comparar la concentración de mercurio en las distintas estaciones en las que hemos capturado los peces, podemos usar:

ggplot(mercurio) +
  geom_boxplot(aes(x = ESTACION, y = CONC, fill = RIO)) +
  labs(title = "Concentración Hg en los peces capturados en varias estaciones",
       x = "Estación",
       y = "Concentración Hg (ppm)")

La última capa modifica algunas anotaciones del gráfico (el título y las etiquetas de los ejes).

Un peligro de los diagramas de cajas es que pueden dar una impresión equivocada si realmente están basados en pocos datos. Vamos a añadir al gráfico anterior los puntos utilizados para calcular las cajas:

ggplot(mercurio, aes(x = ESTACION, y = CONC)) +
  geom_boxplot(aes(fill = RIO)) +
  geom_point() +
  labs(title = "Concentración Hg en los peces capturados en varias estaciones",
       x = "Estación",
       y = "Concentración Hg (ppm)")

Vemos que la caja que parecía mostrar una mayor dispersión (estación 6) corresponde a una estación en la que solo se capturaron cuatro peces.

El gráfico de violín combina un diagrama de cajas con un estimador núcleo (representado dos veces de forma simétrica) de manera que se pueda visualizar simultáneamente la forma de la distribución:

ggplot(mercurio) +
  geom_violin(aes(x = RIO, y = CONC),
              fill = "olivedrab4",
              draw_quantiles = c(0.25, 0.5, 0.75)) +
  labs(title = "Concentración Hg en los peces de dos ríos",
       x = "Río",
       y = "Concentración Hg (ppm)")

1.6.6 Coeficiente de asimetría

Para cuantificar la asimetría de la distribución se usa este coeficiente, que se define como

\[\hat{\gamma} = \frac{1}{n-1} \sum_{i=1}^n \left(\frac{x_i - \bar{x}}{s}\right)^3.\]

Si la distribución es simétrica se espera que se compensen los sumandos positivos con los negativos de forma que \(\gamma \approx 0\). Si la distribución es asimétrica a la derecha predominarán los sumandos muy positivos por lo que \(\gamma>0\). Lo contrario sucede si la distribución es asimétrica a la izquierda.

Al depender únicamente de los datos estandarizados, el coeficiente de asimetría no depende de las unidades de medida.

Para los datos de mercurio en los peces obtenemos los valores siguientes (previamente escribimos una función muy sencilla que lo calcula):

asimetria <- function(x){
  n <- length(x)
  sum((x - mean(x))^3) / ((n-1) * sd(x)^3)
}

mercurio %>% 
  group_by(RIO) %>% 
  summarise(asim_hg = asimetria(CONC),
            asim_longitud = asimetria(LONG),
            asim_peso = asimetria(PESO))
#> # A tibble: 2 x 4
#>   RIO     asim_hg asim_longitud asim_peso
#>   <fct>     <dbl>         <dbl>     <dbl>
#> 1 Lumber    1.48          0.559      1.50
#> 2 Wacamaw   0.939         0.363      1.21

1.7 Relación entre dos variables

En esta sección se introducen las técnicas básicas para describir la relación entre dos variables. Comenzamos por gráficos y después definiremos algunas medidas numéricas para cuantificar lo que se observa en los gráficos.

1.7.1 Diagramas de dispersión

La manera más sencilla de visualizar la relación entre dos variables \(x\) e \(y\) es representar cada observación \((x_i, y_i)\) como un punto. La muestra \((x_i,y_i)\), con \(i=1,\ldots,n\), se corresponde con una nube de puntos que se suele llamar también diagrama de dispersión (o scatter plot en inglés).

ggplot(mercurio) +
  geom_point(aes(x = LONG, y = CONC))

Para diferenciar las observaciones de cada uno de los ríos podemos usar distinto color (o un símbolo distinto si no queremos usar colores):

ggplot(mercurio) +
  geom_point(aes(x = LONG, y = CONC, col = RIO))

ggplot(mercurio) +
  geom_point(aes(x = LONG, y = CONC, shape = RIO))

Para un número moderado de variables puede ser útil dibujar los diagramas de dispersión para cada par de variables

# requiere la extensión GGally

ggpairs(mercurio, columns = 3:5, aes(col = RIO))

Al mirar un diagrama de dispersión conviene fijarse en los siguientes aspectos:

• ¿Hay asociación entre \(X\) e \(Y\)?
• ¿Es la asociación lineal o no lineal? ¿Es positiva o negativa?
• ¿Hay heterocedasticidad? Esto significa que la varianza de la \(Y\) depende de los valores de \(X\). Lo contrario de la heterocedasticidad es la homocedasticidad. Muchos modelos estadísticos asumen que los datos son homocedásticos, por lo que si esta hipótesis no se cumple puede ser necesario aplicar alguna transformación.
• ¿Hay datos atípicos? En un diagrama de dispersión hay atípicos que no lo son respecto a ninguna de las variables por separado, sino porque la relación entre las variables en esos puntos es diferente a la de la mayoría de los puntos.
• ¿Hay grupos homogéneos (clústeres)? La presencia de estos grupos puede indicar la influencia de una variable no observada, de forma análoga a la presencia de varias modas en un histograma.

En los gráficos anteriores se observa una relación positiva entre las variables, que es lineal entre la concentración de mercurio y la longitud o el peso, pero no es lineal entre el peso y la longitud. Además, la relación entre el peso (o la longitud) y la concentración es heterocedástica: para peces grandes la variabilidad en las cantidades de mercurio registradas es mayor que para peces pequeños. En el gráfico entre el peso y la longitud se observa un grupo de peces del río Lumber cuya longitud es mayor que la del resto de acuerdo con su peso. ¿Podrían ser ejemplares de otra especie? En el gráfico de la concentración y el peso también se observan algunos posibles atípicos en el río Lumber para los que la concentración resulta anormalmente baja para el peso que tienen.

1.7.2 Covarianza y correlación

Los gráficos anteriores sirven para determinar visualmente el tipo de relación entre las variables. Para cuantificar el grado de relación lineal se usan dos medidas fundamentales en estadística, la covarianza y la correlación.

Supongamos que tenemos una muestra de datos bidimensionales \((x_1,y_1),\ldots,(x_n,y_n)\) de dos variables \(x\) e \(y\). Se define la covarianza entre \(x\) e \(y\) como \[s_{x,y} = \frac{1}{n-1} \sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y}).\]

La figura 1.5 puede ayudar a interpretar intuitivamente la medida anterior. Se muestran situaciones en las que podemos esperar covarianza positiva, negativa o aproximadamente igual a cero. El vector de medias divide el plano en cuatro cuadrantes. Los puntos situados en el primero y el tercero aportan cantidades positivas a la covarianza mientras que los situados en el segundo y el cuarto aportan cantidades negativas.

Figura 1.5: Interpretación del signo de la covarianza

Otras propiedades elementales de la covarianza son las siguientes:

• La covarianza es simétrica: \(s_{x,y} = s_{y,x}\).
• La covarianza de una variable consigo misma es la varianza de esa variable: \(s_{x,x} = s^2_x\).
• Si definimos los vectores \(x = (x_1-\bar{x},\ldots,x_n-\bar{x})\) e \(y = (y_1-\bar{y},\ldots,y_n-\bar{y})\), entonces \((n-1)s_{x,y} = \langle x, y\rangle\), donde \(\langle \cdot, \cdot\rangle\) denota el producto escalar euclídeo. (Y de igual forma \((n-1)s^2_x = \|x\|^2\), donde \(\|\cdot\|\) es la norma euclídea.)
• Una fórmula alternativa de la covarianza es: \[s_{x,y} = \frac{n}{n-1} \left[\frac{\sum_{i=1}^n x_iy_i}{n} - \bar{x}\bar{y}\right].\]
• La covarianza depende de las unidades en las que medimos las variables. De hecho, si \(\lambda\in\mathbb{R}\), ¿cuánto vale \(s_{\lambda x,y}\)?

Esta última propiedad dificulta mucho la interpretación de la covarianza. De hecho, lo único que podemos interpretar de ella es el signo. Para corregir este problema se normaliza la covarianza dividiendo por el producto de desviaciones típicas. Esto da lugar al coeficiente de correlación entre \(x\) e \(y\): \[r_{x,y} = \frac{s_{x,y}}{s_xs_y}.\]

• Si \(\lambda\in\mathbb{R}\), ¿cuanto vale \(r_{\lambda x,y}\)? ¿Cuánto vale \(r_{x,x}\)?
• En función de normas y productos escalares, si definimos los vectores \(x = (x_1-\bar{x},\ldots,x_n-\bar{x})\) e \(y = (y_1-\bar{y},\ldots,y_n-\bar{y})\), entonces \[r_{x,y} = \frac{\langle x, y\rangle}{\|x\|\cdot \|y\|}.\] Así que podemos interpretar la correlación como el coseno del ángulo que forman \(x\) e \(y\).
• Aplicando la desigualdad de Cauchy-Schwarz, \(|r_{x,y}|\leq 1\). Además, \(|r_{x,y}|=1\) si y solo si existen \(a,b\in\mathbb{R}\) tales que \(y_i = ax_i + b\), para todo \(i = 1,\ldots,n\).

Para calcular la covarianza y la correlación con R basta usar los comandos cov y cor de la siguiente forma:

set.seed(123)
n = 100
x <- rnorm(n)
y1 <- rnorm(n, sd = 10)
cov(x, y1)
#> [1] -0.4372107
cor(x, y1)
#> [1] -0.04953215

# cov y cor también se aplican a matrices
y2 <- 0.8*x + rnorm(n, sd = 0.8)
datos <- cbind(x, y1, y2)

cov(datos)
#>             x         y1         y2
#> x   0.8332328 -0.4372107  0.5769831
#> y1 -0.4372107 93.5063098 -0.1250689
#> y2  0.5769831 -0.1250689  0.9673568
cor(datos)
#>              x          y1         y2
#> x   1.00000000 -0.04953215  0.6426678
#> y1 -0.04953215  1.00000000 -0.0131503
#> y2  0.64266784 -0.01315030  1.0000000

Para conjuntos de datos se pueden usar estos comandos como parte de otros comandos de dplyr, de manera similar a como hemos calculado otras medidas numéricas. Por ejemplo, el siguiente código calcula las correlaciones entre la concentración de mercurio, la longitud y el peso, por separado para cada los peces de cada río:

mercurio %>% 
  group_by(RIO) %>% 
  summarise(cor_hg_long = cor(CONC, LONG),
            cor_hg_peso = cor(CONC, PESO),
            cor_long_peso = cor(LONG, PESO))
#> # A tibble: 2 x 4
#>   RIO     cor_hg_long cor_hg_peso cor_long_peso
#>   <fct>         <dbl>       <dbl>         <dbl>
#> 1 Lumber        0.554       0.381         0.872
#> 2 Wacamaw       0.707       0.704         0.939

1.7.3 Regresión lineal

A veces no solo es necesario cuantificar el grado de relación entre dos variables sino aproximar una de las variables en función de la otra. En ese caso tenemos que determinar lo que se conoce como recta de regresión. El objetivo en este apartado es analizar la relación existente entre dos variables, \(x\) e \(y\), de forma que podamos predecir o aproximar el valor de la variable \(y\) a partir del valor de la variable \(x\). A la variable \(y\) se llama variable respuesta. A La variable \(x\) se llama variable regresora o explicativa.

En un problema de regresión (a diferencia de cuando calculamos el coeficiente de correlación) el papel de las dos variables no es simétrico.

Recta de mínimos cuadrados

Una recta \(y = \beta_0+\beta_1 x\) adecuada es aquella para la que las diferencias \(y_i - (\beta_0+\beta_1 x_i)\) para \(i=1,\ldots,n\) tienden a ser pequeñas (sin tener en cuenta el signo).

La recta de regresión de mínimos cuadrados viene dada por los valores \(\hat{\beta}_0\) y \(\hat{\beta}_1\) para los que se minimiza la función \[\phi(\beta_0,\beta_1)=\sum_{i=1}^n [y_i - (\beta_0+\beta_1 x_i)]^2\] Dado que esta función es convexa, para encontrar la recta basta determinar los valores para los que se anula su gradiente, es decir, resolver el siguiente sistema de ecuaciones (la primera ecuación resulta de derivar respecto a \(\beta_0\) y la segunda respecto a \(\beta_1\)):

\[\begin{align*} \sum_{i=1}^n [y_i - (\hat{\beta}_0+\hat{\beta}_1 x_i)] &= 0 \\ \sum_{i=1}^n x_i[y_i - (\hat{\beta}_0+\hat{\beta}_1 x_i)] &= 0. \end{align*}\]

De la primera ecuación se deduce fácilmente que la recta que buscamos pasa por el punto \((\bar{x},\bar{y})\). Un ejercicio sencillo permite calcular también la pendiente, de manera que la recta de mínimos cuadrados es: \[y - \bar{y} = \hat{\beta}_1 (x -\bar{x}),\] donde la pendiente viene dada por \[\hat{\beta}_1 = \frac{s_{x,y}}{s^2_x} = r_{x,y} \frac{s_y}{s_x} = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}.\]

Residuos y valores ajustados

El valor pronosticado o previsto (también llamado valor ajustado) de \(y_i\) dado \(x = x_i\) es \(\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\).

Las diferencias entre los valores reales de la variable respuesta y los valores ajustados, \(e_i = y_i - \hat{y}_i\), se denominan residuos. La recta de mínimos cuadrados es aquella para la que el valor de \(\sum_{i=1}^n e^2_i\) es el mínimo posible. Además, las dos ecuaciones que hemos tenido que resolver para calcularla permiten deducir que \(\sum_{i=1}^n e_i = 0\) y que \(\sum_{i=1}^n x_i e_i = 0\). Equivalentemente, los residuos tienen media cero y están incorrelados con la variable regresora (¿por qué esto es equivalente a lo anterior?).

Una manera de medir el error cometido en la aproximación dada por la recta de regresión es mediante la llamada varianza residual \[s^2_R = \frac{\sum_{i=1}^n e_i^2}{n-2}.\] Podemos interpretar esta cantidad como la incertidumbre que aún se tiene sobre la variable respuesta una vez observada la variable regresora. Puede demostrarse (ejercicio) que \[\frac{\sum_{i=1}^n e_i^2}{n-1} = s^2_y (1-r^2_{x,y}),\] y por tanto \[\frac{n-2}{n-1} s^2_R = s_y^2 (1-r_{x,y}^2),\] y se verifica que \[r_{x,y}^2 \approx 1 - s^2_R/s_y^2.\] Esta última propiedad permite interpretar el coeficiente de correlación al cuadrado, \(r_{x,y}^2\), como la fracción de la varianza de \(y\) que ya queda explicada después de conocer \(x\). En este contexto, al valor del coeficiente de correlación al cuadrado también se le llama coeficiente de determinación.

Cálculo y representación de la recta

El comando básico para calcular la recta de mínimos cuadrados es lm. Se usa de la siguiente forma:

# Cálculo de la ecuación de la recta de mínimos cuadrados
lm(CONC ~ LONG, data = mercurio)
#> 
#> Call:
#> lm(formula = CONC ~ LONG, data = mercurio)
#> 
#> Coefficients:
#> (Intercept)         LONG  
#>    -1.13165      0.05813

Así obtenemos los valores de \(\hat\beta_0\) (intercept) y \(\hat\beta_1\), usando los peces de los dos ríos conjuntamente. Para obtener las ecuaciones de la recta por separado podemos filtrar previamente los datos, tal y como hemos visto en una sección anterior:

mercurio_lumber <- mercurio %>% 
  filter(RIO == 'Lumber')

lm(CONC ~ LONG, data = mercurio_lumber)
#> 
#> Call:
#> lm(formula = CONC ~ LONG, data = mercurio_lumber)
#> 
#> Coefficients:
#> (Intercept)         LONG  
#>    -0.62387      0.04318

La representación gráfica se puede obtener mediante geom_smooth(). Se usa el argumento method = 'lm' para que ajuste una recta en lugar de una función no lineal y el argumento se= FALSE para que no represente gráficamente las bandas de confianza:

# Representación gráfica
ggplot(mercurio, aes(x = LONG, y = CONC, col = RIO)) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE) +
  scale_colour_brewer(palette="Set1")

Hemos usado la última línea para cambiar la paleta de colores por defecto. Para más información sobre cómo seleccionar los colores en ggplot2 se puede consultar esta página.

Otra posibilidad es representar un gráfico por separado para cada río. Esto se consigue fácilmente añadiendo una capa facet (que es la jerga que usa ggplot2 cuando se trata de hacer gráficos separados para cada uno de los niveles de un factor):

ggplot(mercurio, aes(x = LONG, y = CONC)) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE) +
  facet_wrap(~RIO)

Si queremos que los gráficos aparezcan en columna usamos el argumento ncol:

ggplot(mercurio, aes(x = LONG, y = CONC)) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE) +
  facet_wrap(~RIO, ncol = 1)

1.8 Ejercicios

Ejercicio 1.1 Realiza un análisis descriptivo de los datos británicos de ingresos familiares contenidos en el fichero Datos-ingresos.txt. En concreto, calcula los estadísticos de posición o tendencia central, las medidas de dispersión, representa un diagrama de cajas y un estimador del núcleo de la función de densidad. Comenta los resultados.

Ejercicio 1.2 Representa en el mismo gráfico los diagramas de cajas correspondientes a la variable Largo del fichero tortugas.txt para los ejemplares hembra y para los ejemplares macho. Emplea colores distintos para los dos diagramas.

Ejercicio 1.3 Los datos del fichero kevlar.txt corresponden al tiempo hasta el fallo (en horas) de 101 barras de un material utilizado en los transbordadores espaciales, llamado Kevlar49/epoxy, sometidas a un cierto nivel de esfuerzo. Los datos han sido tomados de Barlow et al. (1984).

• Calcula las principales medidas descriptivas numéricas de estos datos.
• Representa un diagrama de cajas.
• Representa un histograma con un número de clases apropiado.
• Estudia la presencia de datos atípicos en la muestra. Si hay datos atípicos, suprímelos y repite todos los apartados anteriores. Compara los resultados obtenidos.

Ejercicio 1.4 El paquete gapminder contiene un fichero de datos de población, esperanza de vida y renta per cápita de los países del mundo entre 1952 y 2007. Instala el paquete y lleva a cabo los siguientes gráficos:

1. Un histograma de la esperanza de vida en 2007 de los países de Europa.
2. Diagramas de cajas con las esperanzas de vida de cada continente en el año 1952.
3. Un diagrama de dispersión de la renta per cápita y la esperanza de vida de cada país en el año 2007.
4. Mejora el gráfico anterior representando cada punto con un color diferente en función del continente al que pertenece cada país y representando la renta per cápita en una escala logarítmica (Indicación: mira esta documentación)

Ejercicio 1.5 Calcula el diagrama de dispersión de las dos variables correspondientes al peso y a la circunferencia de abdomen que aparecen en el fichero Datos-bodyfat.txt. Calcula la recta de regresión y el coeficiente de correlación. Comenta los resultados.

Ejercicio 1.6 El fichero star.txt contiene la temperatura y la intensidad de la luz en un conjunto de estrellas. Calcula y representa la recta de mínimos cuadrados para explicar la temperatura en función de la intensidad de la luz. Comenta el resultado.

1.9 Referencias

Los capítulos 15 y 16 de Dekking et al (2005) contienen una breve introducción a la descripción de datos, que incluye también los estimadores de núcleo de la densidad, lo que no es habitual en este tipo de libros. Grolemund y Wickham (2017) es una buena introducción al tratamiento de datos básico con R (usando tidyverse). Otro libro en la misma línea es Irizarry (2020). Wilke (2019) contiene muchas ideas razonables acerca de cómo visualizar conjuntos de datos, sin que el libro esté ligado a ningún software en concreto. Por el contrario, Healy (2018) aborda el problema a través del uso de ggplot2 como herramienta de visualización.

• Dekking, F. M., Kraaikamp, C., Lopuhaä, H. P. y Meester, L. E. (2005). A Modern Introduction to Probability and Statistics: Understanding why and how. Springer.
• Grolemund, G. y Wickham, H. (2017). R for Data Science. O’Reilly | Versión en castellano
• Healy, K. (2018). Data Visualization, a practical introduction. Princeton University Press.
• Irizarry, R.A. (2020). Introduction to Data Science. CRC Press.
• Wilke, C. O. (2019). Fundamentals of data visualization. O’Reilly