Ciencia de datos con R

Table of Contents

1 Fundamentos de R

Estos son los apuntes del curso Data Science with R ofrecido por la universidad de Harvard a través de la plataforma HarvardX.

1.1 Lo más básico

1.1.1 Instalación de R

Instalar R en Debian es sencillo: apt install r-base r-recommended.

Se inicia un terminal de R con la orden R.

1.1.2 Instalar y cargar paquetes

Instalar paquetes aumenta la funcionalidad

install.packages("nombre_del_paquete")

Se pueden instalar en grupos, más de uno a la vez usando un vector como entrada de la orden de instalación, así:

install.packages(c("paquete1", "paquete2", "paquete3"))

En Debian, R intenta instalar los paquetes en /usr/lib y ese directorio no es escribible por el usuario. Se puede definir una ruta alternativa (p. ej. en /home/usuario/.R/library) pero para que R pueda compilar e instalar el código de los paquetes deseados es necesario instalar en el sistema el paquete libssl-dev.

Es interesante guardar un script con la instalación de todos los paquetes usados normalmente, de manera que es sencillo reinstalarlos cada vez que se quiera utilizar una instalación base de R o otro ordenador en el que no se haya trabajado anteriormente. La lista de paquetes instalados se puede recuperar con ìnstalled.packages().

Solo es necesario instalar los mismos paquetes otra vez si se utiliza R con otro usuario en el mismo ordenador u otro ordenador distinto. O, por supuesto, si se borra el directorio que contiene los paquetes descargados y compilados.

También es bueno empezar un script con la instalación y carga de los paquetes necesarios para su ejecución.

Una vez instalado, un paquete que no pertenezca a la instalación base de R necesita ser cargado con la instrucción library.

library(nombre_del_paquete)

La instalación se hace una vez, pero la carga del paquete se hace siempre que haga falta (una vez por sesión o script). Nótese la falta de comillas dentro de los paréntesis una vez que el paquete está instalado. Cuando se intenta cargar un paquete y se muestra un error, lo más probable es que sea necesario instalarlo primero.

1.1.3 Objetos

Los objetos son cosas guardadas en R. Se asignan valores a variables y eso es un objeto:

a <- 1

En este caso, a toma el valor 1 y se puede operar con a igual que con un 1: a + 1 = 2.= Pero los objetos pueden ser mucho más complicados, como funciones. Los valores se asignan con la flecha <-, aunque también sirve el signo igual =. No se recomienda usar = para evitar confusiones, solo la flecha.

Escribir tan solo el nombre de la variable hace que R muestre su valor, por ejemplo:

a
[1] 1

muestra un 1 como valor de la variable y, entre corchetes, que este es el primero de los datos mostrados. A veces, un vector consta de muchos datos en una lista y, al principio de cada línea de esta lista se muestra entre corchetes el número de orden del primer elemento de esa línea.

Si la variable (el objeto) no tiene un valor asignado, R mostrará un mensaje de error. Es decir, no se pueden recuperar valores de objetos aún no definidos. Por ejemplo:

a <- 1
x
Error: object 'x' not found.

Para definir una variable hay que seguir unas reglas sencillas: 1. Han de comenzar con una letra. 2. No pueden contener espacios. 3. No pueden coincidir con variables predefinidas en R (como pi o install.packages).

Una buena costumbre es usar nombres que sean explicativos y sustituir los espacios con guiones bajos (primera_solucion, datos2, etc.).

La orden ls() elabora una lista con los objetos que R tiene guardados en el espacio de trabajo (ver más abajo). Cada vez que se define una variable o se asignan valores a un vector, estos elementos aparecen en la salida del comando ls().

1.1.4 Espacios de trabajo

Cada vez que se define una variable u otro objeto, se está modificando el espacio de trabajo (workspace). Este puede ser guardado al finalizar la sesión y cargado al comienzo de la siguiente (o deshechado). La orden ls() muestra una lista de los objetos guardados en el espacio de trabajo. Hay órdenes (como save o load) para guardar y cargar espacios de trabajo, pero aún no sé cómo funcionan, no parece fácil al principio.

Estos objetos se muestran en RStudio (in IDE para R) en la pestaña Environment.

1.2 Funciones

Una vez definidas variables, el análisis de los datos se puede definir como una serie de funciones aplicadas a los datos. Algunas de estas funciones ya están predefinidas en R.

Las órdenes usadas antes como install.packages, library son funciones. Se pueden cargar más funciones que las predefinidas instalando y cargando paquetes.

Cuando se evalúa una función, que es lo mismo (más o menos) que ejecutarla, hay que incluir paréntesis aunque no haya nada entre ellos. Si se necesita pasarle parámetros a la función, han de incluirse entre los paréntesis. Si no se usan los paréntesis, R muestra el código fuente de la función.

Las funciones incorporan un documento de ayuda para aprender lo que hacen y cómo: qué parámetros son necesarios, etc. Por ejemplo, help("log") muestra la ayuda de la función que calcula el logaritmo de un número: log(x). Se consigue el mismo resultado con el signo de interrogación: ?log.

En la ayuda se especifican los parámetros necesarios, obligatorios u opcionales, y el orden en el que deben aparecer. Por eso, si se respeta el orden no es necesario definirlos antes. Por ejemplo:

log(x = 8, base = 2)

es equivalente a

log(8,2)

porque los parámetros están en el orden esperado. Pero también funciona lo siguiente:

log(base = 2, x = 8)

donde los parámetros no están en orden pero definidos.

La función args("nombre_de_la_funcion") devuelve los parámetros de una función; sirve para tener a mano un recordat orio rápido sin tener que acceder al documento completo de la ayuda.

Hay otros objetos predefinidos, como constantes matemáticas (pi (pi), infinito (inf) u otras) o incluso conjuntos de datos completos. Estos conjuntos de datos de ejemplo se pueden consultar con la instrucción data().

Las funciones pueden recibir tanto datos (numéricos o no, depende de la función) como variables, que serán interpretadas como los datos asignados a ellas. Por eso, si se asigna x <- 8, tanto log(8) como log(x) darán el mismo resultado.

Las funciones se pueden anidar, y entonces se usa el resultado de unas como parámetro para las siguientes. El orden es siempre de dentro afuera: se evalúan primero las más interiores, las más profundamente anidadas, y luego las más exteriores, las que contienen a las demás. log(sqrt(9)) calcula primero la raíz de 9 y luego el logaritmo natural de ese resultado. Ya se usó la anidación de funciones en un apartado anterior al instalar paquetes de R que venían dados por una lista asignada a un vector.

1.2.1 Escribir programas

Aunque definir variables y escribir fórmulas y funciones para operar con ellas desde luego que funciona, es mucho más eficaz escribir todo en un documento de texto y luego ejecutarlo. De esta manera, tan solo cambiando el valor de las variables pueden obtenerse diferentes resultados.

En un programa se pueden escribir comentarios, es decir, líneas que no serán ejecutadas y sirven para documentar el código, como por ejemplo:

# Asignar valores a las variables
a <- 1
b <- 2
c <- 3

# Fórmula para ser ejecutada
resultado <- sqrt(a+b+c)

# Mostrar el resultado:
resultado

Para escribir un programa o script, simplemente se escribe su contenido en un archivo de texto y se guarda con extensión R: archivo.R1. Luego se ejecuta en una sesíón de R con la orden source(ruta/al/archivo.R).

Es muy recomendable escribir los programas en archivos de texto ejecutables porque permite revisar el trabajo, repetirlo cambiando parámetros, corregir errores fácilmente, etc. También es más sencillo organizar los programas y compartirlos.

Las líneas que comienzan con almohadillas (#) no son leídas por el programa y por lo tanto no se ejecutan. Por lo tanto valen para insertar comentarios como se mencionó antes y hacer el código más fácil de entender. Es algo muy recomendable.

1.3 Tipos de datos

La función class devuelve el tipo de objeto que se le pasa como parámetro:

a <- 1
class(a)
[1] numeric

Aunque se pueden asociar variables a números, no es la mejor ni más eficiente manera de manejar datos, lo más habitual es organizar un conjunto de datos (data frame). Así se pueden manejar diferentes tipos de datos en un solo objeto.

En el ejemplo se usa el conjunto de datos proporcionado para analizar el número de asesinatos en EE.UU. Es necesario tener acceso a él, cargarlo y comprobar de qué tipo de objeto se trata.

install.package("dslabs")
library(dslabs)
data("murders")
class(murders)

[1] "data.frame"

Para saber más acerca del conjunto de datos se usa str (de structure: str(murders)). Esta función devuelve una vista de la estructura de los datos con el número de observaciones (filas, entrevistados, etc.), el número de variables (columnas, preguntas, opciones, etc.) y los nombres de las variables con los primeros valores de cada una a modo de ejemplo.

La función head (head(murders)) muestra los primeros 6 elementos de cada variable en forma de tabla, con las variables como columnas y las observaciones (las seis primeras( en las filas.

1.3.1 Acceder a las variables

Para poder saber algo acerca de los datos contenidos en el objeto hay que estudiar las variables por separado y relacionar unas con otras. Para acceder a los datos de las variables se usa el signo $ en el nombre del conjunto de datos: murders$population muestra los datos contenidos en esa variable, 51 números en total, uno por estado (los estados son, en este caso, las observaciones). Los datos mostrados (y esto es importante) siguen el orden de las observaciones tal cual están registradas en el conjunto de datos, fila por fila.

Los nombres de las variables aparecen como una lista en la salida de la función str() y como las cabeceras de las columnas en la salida de la función head(), pero también pueden consultarse específicamente con la función names().

Pero, empezando por el principio, lo más sencillo es entender primero los vectores.

1.3.2 Vectores: numéricos, alfanuméricos y lógicos

Tipos de vectores (los básicos y más usuales): 1. Numéricos. 2. Alfanuméricos o texto. 3. Lógicos. 4. Factores.

Ciertos objetos, como la información contenida en la variable murders$population, no constan de un valor numérico sino de muchos (51 en el ejemplo). Este tipo de objetos es un vector. Teóricamente, un número es un vector de longitud 1.

La función length informa de la cantidad de datos que hay en un vector. En el siguiente ejemplo se asigna el vector a un objeto y se comprueba la longitud de ese vector. También se puede comprobar directamente: ambas instrucciones son equivalentes.

pop <- murders$population
length(pop)
length(murders$population)
[1] 51

Los vectores pueden almacenar números, como murders$population. La función class aplicada a este vector devolverá numeric. Estos son vectores numéricos.

También pueden almacenar texto, como el vector formado por la columna (o variable) state en este mismo conjunto de datos. Así, class(murders$state) informará de que el tipo de vector es alfanumérico, o sea, texto: character.

Un tercer tipo de vectores es el vector lógico, que puede contener solo los valores TRUE o FALSE.

z <- 3 == 2
z
 [1] FALSE

class(z)
 [1] "logical"

En este ejemplo se ve el operador == que «pregunta» si 3 es igual a 2. Usando = se asigna un valor a una variable.

Los factores son vectores alfanuméricos pero que almacenan solo un número limitado de cadenas diferentes. En el conjunto de datos del ejemplo, el vector murders$regions es un vector de factores porque almacena categorías, no diferentes datos (solo cuatro regiones diferentes en todo el conjunto de datos guardadas como 2 4 4 2 4 4 1 2… según se ve en str(murders)).

La función levels muestra las categorías de que se compone el vector de factores.

levels(murders$region)

Este ejemplo devuelve cuatro valores: Northeast, South, North Central y West.2 No importa el número de filas (observaciones) ni el número de veces que aparezca cada una de las regiones en el conjunto de datos, estos son todos los tipos que existen. Anidando funciones, se puede saber el número de niveles (valores) únicos de la variable de tipo factor con la función length() así: length(levels(murders$regions)). Esta función anidada devuelve el valor 4: hay cuatro diferentes factores en todo el conjunto de datos.

Las confusiones entre factores y vectores alfanuméricos son una fuente de problemas en los programas.

1.3.3 Trabajar con un conjunto de datos (data frame)

Un data frame o conjunto de datos es, básicamente, una tabla en la que pueden convivir datos numéricos, alfanuméricos, lógicos y de cualquier tipo.

Para los ejemplos a continuación se usará un conjunto de datos de ejemplo que consiste en los ratios de asesinatos por arma de fuego en los EE.UU. para cada estado, junto a su población y otros detalles. Para acceder a este conjunto de datos hay que instalar el paquete correspondiente y cargarlo para que sus datos sean accesibles. Luego cargar el conjunto de datos en sí.

install.packages("dslabs")
library(dslabs)
data(murders)

Una vez que el paquete ha sido instalado no es necesario volver a hacerlo en el mismo ordenador. A partir de ese momento es suficiente con cargar la librería y el conjunto de datos en cada sesión.

La función str(objeto) es útil para analizar la estructura del objeto, pues muestra un resumen. Por ejemplo:

str(murders)

mostrará que es un conjunto de datos (data.frame), el número de observaciones (filas) y variables (columnas), y un resumen del contenido de estas variables.

La función names devuelve los nombres de las columnas de un conjunto de datos. Serían los nombres de las variables o columnas contenidas en cada observación o fila. En el caso del ejemplo, sería names(murders).

Más información se puede extraer con head(murders). Esta función muestra las seis primeras filas del conjunto de datos.

Cada variable tiene una clase, que puede ser numérica, alfanumérica o lógica. Se puede consultar la clase con la función class(). Pero para extraer los datos de una sola variable, es necesario usar el símbolo $ (llamado accesor 'accessor') antes del nombre de la variable. Así class(murders$abb) muestra la clase de datos en el vector formado por los valores de la variable abb dentro del conjunto de datos murders. Como se dijo, para conocer los nombres de las variables incluidas en un conjunto de datos (las columnas), se puede usar la función names.

Los valores de una variable en un conjunto de datos se pueden extraer también con el uso de corchetes: x <- murders[["population"]].

a <- murders$population
b <- murders[["population"]]
identical(a, b)

Si en lugar de corchetes dobles se usan corchetes sencillos, el objeto al que se asocie será un subconjunto de datos de clase data.frame en lugar de un vector.

Dentro de un conjunto de datos puede haber variables de tipo factor. Son valores numéricos de un rango determinado que se asocian con un nombre o valor, como el valor "Region" en el ejemplo (solo cuatro regiones diferentes en todo el conjunto de datos guardadas como 2 4 4 2 4 4 1 2… según se ve en str(murders)). Los valores (niveles) que puede tomar la variable de tipo factor se pueden consultar con la función levels().

levels(murders$region)

Este ejemplo devuelve cuatro valores: Northeast, South, North Central y West. No importa el número de filas (observaciones) ni el número de veces que aparezca cada una de las regiones en el conjunto de datos, estos son todos los tipos que existen. Anidando funciones, se puede saber el número de niveles (valores) únicos de la variable de tipo factor con la función length() así: length(levels(murders$regions)). Esta función anidada devuelve el valor 4.

¿Cuántos estados hay por región? La función table() muestra el número de observaciones en cada nivel de la variable de tipo factor, con su nombre asociado. En este caso, table(murders$region).

1.3.4 Vectores

Un vector es la unidad básica de almacenamiento de datos en R. Es una list de datos. Los conjuntos de datos complejos se suelen poder separar en vectores individuales. Por ejemplo, en el conjunto murders, cada una de las variables es a su vez un vector (cada uno de una clase diferente, etc.).

La forma más sencilla de crear un vector es con la función c, de concatenate. Se pueden crear vectores de clase numérica tanto como alfanumérica:

codes <- c(380, 124, 818)
country <- c("italy", "canada", "egypt")

Es importante saber que se pueden otorgar nombres a los elementos de un vector numérico, y se puede hacer de varias maneras: definiéndolo en el mismo vector, o creando objetos con datos y objetos con nombres para los datos, y uniéndolos después. Ejemplos:

codes <- c(italy = 380, canada = 124, egypt = 818)
codes <- c("italy" = 380, "canada" = 124, "egypt" = 818)

codes <- c(380, 124, 818)
country <- c("italy","canada","egypt")
names(codes) <- country

Otra manera de crear un vector es con la función seq(a, b) (viene de sequence 'secuencia'). Esta función muestra los enteros consecutivos de a a b. Aunque el comportamiento por defecto es mostrar los enteros de uno en uno, un tercer parámetro permite establecer el valor de la diferencia entre elementos. Así, seq(1, 10, 2) mostrará los impares desde 1 (el comienzo) hasta 10 (el final): 1, 3, 5, 7, 9.

También es válido el siguiente atajo: [1:10] muestra los números enteros consecutivos entre 1 y 10.

1.3.5 Acceder a elementos de un vector

Los corchetes sirven para acceder a elementos específicos dentro de un vector. Por ejemplo, para los vectores anteriores, los siguientes códigos muestran los datos, dentro del vector codes definido antes, de la posición 2, de las posiciones 1 y 3, de las posiciones 1 a 2 (no 1 Y 2), los correspondientes al nombre canada y a los nombres egypt e italy3:

codes[2]
codes[c(1,3)]
codes[1:2]
codes["canada"]
codes[c("egypt","italy")]

Nótese que el parámetro que se le pasa al nombre del objeto (codes) es en sí mismo un vector también, por lo que han de «aglutinarse» con la función concatenate: c(1. 3).

1.4 Coerción (forzado) de vectores

La coerción es el intento de R de ser flexible con los tipos de datos utilizados. Cuando los datos introducidos no concuerdan con lo esperado se intenta entender lo que se quiso decir antes de mostrar un mensaje de error.

Coerción4 significa lo mismo que 'extorsión' o 'coacción': Presión ejercida sobre alguien para forzar su voluntad o su conducta; y también Represión, inhibición, restricción. Coerción en este contexto se refiere a la manera en que R restringe el uso de un solo tipo de datos en un vector. Si hay datos numéricos, son numéricos; si hay datos alfanuméricos, son alfanuméricos; si son lógicos, son lógicos. Pero cuando hay varios tipos mezclados la norma es:

  • Si hay cadenas alfanuméricas, R lo convertirá todo a cadenas alfanuméricas.
  • Si hay valores lógicos mezclados con valores numéricos, R convertirá los valores lógicos a numéricos asimilando TRUE a 1 y FALSE a 0.

Por eso se pueden hacer operaciones con vectores lógicos, porque los TRUE son tratados como 1 y los FALSE tratados como 0. De esta manera, se pueden conocer la suma de TRUE en un vector o la media de los valores (1 sería la media si todos los valores fueran TRUE y 0 sería la media si todos los valores fueran FALSE; si fueran mitad y mitad, la media sería 0.5).

De manera que el código x <- c(1, "Canadá", 3), aunque tiene valores numéricos, será un vector alfanumérico, pues tanto 1 como 3 pueden ser tratados como caracteres en lugar de números. En efecto, class(x)=devolverá ="character", y los valores almacenados en x son "1" "Canadá" "3", con los números entre comillas porque son tratados como caracteres.

R tiene funciones para forzar un tipo específico de coerción. La función as.character() fuerza que los números sean considerados caracteres, es decir, cadenas alfanuméricas.

  x <- 1:5
  y <- as.carachter(x)
  y
#  [1] "1" "2" "3" "4" "5"

La operación inversa se logra con la función as.numeric. as.numeric(y) devuelve 1 2 3 4 5.

Esto es muy útil cuando los datos son presentados en formas que hacen que R interprete números como caracteres o al revés (códigos postales, por ejemplo).

1.4.1 Datos que faltan (NA)

En R, los datos que faltan en un vector, algo muy habitual en conjuntos de datos reales, se muestran como NA, de not available. Cuando R intenta forzar un vector pero no puede, aparecen datos de esta manera. Por ejemplo, intentado convertir letras en números:

x <- c("1", "b", "3")
as.numeric(x)
 [1] 1 NA 3
 Warning message:
 NAs introduced
 by coercion

Lo que pasa es que los vectores deben guardar el mismo tipo de datos, como se dijo. Sin embargo, R trata de entender lo que "significan" los datos y ajustarlos en la medida de lo posible. Por ejemplo, convertirá el vector c(1, "canada", 3) en un vector alfanumérico "1" "canada" "3". Esto en ocasiones puede causar confusión cuando se cometen errores y los datos son convertidos por fuerza (coerced) en datos de otro tipo.

Las funciones as.character() y as.numeric() convierten vectores a sus respectivos tipos de datos. Cuando la conversión no es posible, como en el ejemplo anterior, en el que "canada" no puede ser convertido a números, R asigna el valor NA (Not Available) al dato que falta y avisa de que la coerción ha introducido valores "NA".

Para operar con el conjunto de datos obviando los valores NA se puede usar el operador na.rm. Por ejemplo:

library(dslabs)
data(na_example)

### Todos los resultados son NA:
mean(na_example)
sd(na_example)

### Se hacen las operaciones normal
### obviando los NA:
mean(na_example, na.rm = TRUE)
sd(na_example, na.rm = TRUE)

1.5 Ordenar datos de un vector

Hay varias funciones útiles para ordenar datos en un vector, o también podría decirse que para ordenar un vector en función de sus datos.

Las funciones sort, order y rank proporcionan utilidades para ordenar los datos dentro de un vector. En realidad no modifican nada, solo muestran los datos con una ordenación determinada.

1.5.1 Función sort

Esta función simplemente ordena los datos de menor a mayor. De manera que si se asocia un vector a un objeto y luego se requiere el primer elemento de ese objeto se obtiene el valor más bajo de todo el vector. Ejemplo:

library(dslabs)
data(murders)
pop <- murders$population
pop <- sort(pop)
pop[1]

1.5.2 Función order

Esta función devuelve el vector de índices necesario para ordenar el vector en la función sort. Es decir, es un vector compuesto por los números de orden necesarios para ordenar los elementos del vector sort. Un ejemplo:

Vector: 23, 1, 56, 57, 12
Sort: 1, 12, 23, 56, 57
Order: 2, 5, 1, 3, 4

Los datos de vector están ordenados así: el segundo, el quinto, el primero, el tercero y finalmente el cuarto (de menor a mayor). Es útil para encontrar los números de fila que cumplan ciertas propiedades, como la fila correspondiente al estado con la menor población:

library(dslabs)
data("murders")
pop <- murders$population
ord <- order(pop)
ord[1]

En realidad esto mismo se puede hacer con la función which.min, que busca la posición (la fila, el índice) del valor mínimo de un vector: which.min(murders$population).

1.5.3 Función rank

La función rank=devuelve el vector formado por las posiciones de los elementos en el vector ordenado por la función =sort. Es decir, es un ranking de cada dato, posición a posición. Con un ejemplo se ve más claro:

data vector sort order rank
31 4 2 3
4 15 3 1
15 31 1 2
92 65 5 5
65 92 4 4

En la columna rank el primer elemento del vector de datos (31) es el tercero más bajo (es el tercero en el vector ordenado sort); el segundo elemento del vector de datos (4) es el menor, el más bajo, es el primero en el vector ordenado sort; etc.

1.6 Operaciones con vectores

Es útil ordenar los datos, extraer el valor máximo o el mínimo, pero para sacar partido a un conjunto de datos, suele ser necesario manipular los vectores en base a otros vectores.

Por ejemplo, en el apartado anterior se vio que California es el estado con mayor número de asesinatos: murders$state[which.max(murders$total)]. Pero también es el estado con mayor población total: `murders$state[which.max(murders$population)] (más de 37 millones de habitantes).

Por lo tanto, para hacernos una idea más realista de lo peligroso (o no) que es el estado de California deberíamos medir el número de asesinatos per cápita.

Las operaciones con vectores (y entre vectores) se hacen elemento por elemento. Es decir, si se multiplica el vector c(1, 3, 5, 7, 9) por 2 se obtiene el vector 2, 6, 10, 14, 18, que es cada uno de los elementos multiplicado por 2. Pero cuando las operaciones se hacen con dos vectores de igual longitud, las operaciones se siguen haciendo elemento por elemento, de manera similar a una hoja de cálculo (la primera celda de la primera columna se suma con la primera celda de la segunda columna, y así en cada fila). Por eso, un código como el siguiente:

library(dslabs)
data("murders")
murder_rate <- murders$total/murders$population*100000

dará como resultado una nueva variable (o columna) murderrate (el ratio de asesinatos por 100,000 habitantes en cada estado) en el conjunto de datos con los valores, para cada estado, de la población total dividida por la población total del estado y multiplicado por 100,000.

1.7 Indexación de vectores

Las posibilidades en R a la hora de manipular vectores son muy potentes y útiles. Por ejemplo, se pueden manipular vectores basándose en características de otros vectores.

Siguiendo con el ejemplo, ahora que hemos calculado el índice de asesinatos por cada 100,000 habitantes para cada estado (murder_rate), lo usaremos para escoger un estado al que mudarnos si, por ejemplo, venimos de Italia, que tiene un índice de asesinatos por cada 100,000 habitantes de 0.71. Para ello, usaremos una característica muy potente de R, que es indexar vectores basándose en valores lógicos (TRUE/FALSE). Esto se aplica a las operaciones lógicas, en las que la comparación con un término aplica la operación a cada uno de los elementos del vector. Está más claro con el ejemplo a continuación. El siguiente código:

index <- murder_rate < 0.71

creará un vector lógico en el que habrá valores TRUE en las observaciones (filas) en las que el índice de asesinatos sea menor que 0.71 y FALSE en las que sea igual o mayor. El operador <= dará el mismo resultado pero será TRUE cuando el valor de «murderrate» sea menor o igual que 0.71. Pero esta es una lista de «trues» y «falses» que no dice mucho. Pero filtrar los nombres basándose en esta lista es mucho más útil: murders$state[index] muestra únicamente los valores Hawaii, Iowa, New Hampshire, North Dakota y Vermont.

El vector lógico index es forzado ('coerced') en uno numérico en el que FALSE es 0 y TRUE es 1. Por eso, sum(index) da como resultado 5, porque suma todo y solo hay 5 unos, los correspondientes a los estados de Hawaii, Iowa, NH, ND y Vermont.

Ahora ya sabemos qué estados tienen un índice de asesinatos menor o igual que el de Italia. Ahora bien, si además queremos mudarnos a un estado en el oeste será necesaria una nueva operación de filtrado.

Digamos que el límite de asesinatos por 100,000 habitantes es 0.71 y que la región es el oeste. Por lo tanto, pueden definirse las siguientes condiciones:

safe <- murder_rate <= 0.71
west <- murders$region == "West"
# El operador "&" creará un vector lógico en el que serán "TRUE" los
# elementos que satisfagan AMBAS condiciones
index <- safe & west

murders$state[index]
[1] Hawaii

1.7.1 Funciones útiles para la indexación

Tres nuevas funciones muy útiles para manipular vectores son which, match e %in%.

  1. Which

    which devuelve las entradas (el número de índice, el número de fila) de un vector lógico que son ciertas (es decir, son TRUE).

    x <- c(FALSE, TRUE, FALSE, TRUE, TRUE, FALSE)
    which(x)
    [1] 2 4 5
    

    Por ejemplo, para ver el índice de asesinatos de Massachusetts, se puede asignar el valor de which a una variable basándose en un operador lógico, así:

    index <- which(murders$state == "Massachusetts")
    index
    [1] 22
    murder_rate[index]
    [1] 1.802
    

    Por supuesto, no es la única manera de lograrlo, con lo que se sabe hasta ahora, se puede conocer un índice de asesinatos en concreto filtrando el vector con el operador lógico directamente: index <- murders$state = "Massachusetts"= y, a continuación el filtrado murder_rate[index]. O, también, directamente: murder_rate[murders$state = "Massachusetts"]=.

  2. Match

    Busca entradas en un vector y devuelve los índices (o filas) en las que se puede acceder a esos datos. Como parámetro se puede pasar, a su vez, un vector con datos que se quieren comprobar. Un ejemplo: para conocer el índice de asesinatos de varios estados en lugar de uno solo se puede crear un vector con «el cruce» de la información de dos vectores: el vector de nombres y el vector de estados a consultar. Así:

    index <- match(c("New York", "Florida", "Texas"), murders$state)
    index
    [1] 33 10 44
    murder_rate[index]
    [1] 2.667960 3.398069 3.201360
    # Comprobación de si los índices dados por la función equivalen a los
    # estados que queremos consultar
    murders$state[index]
    [1] "New York" "Florida" "Texas"
    
  3. %in%

    Esta función lo que hace es comprobar si cada elemento de un vector existe en otro vector dado. Una muestra: si se define x como las seis primeras letras x <- c("a", "b", "c", "d", "e") e y como otra secuencia de letras diferente y <- c("a", "d", "f"), la función %in% devolverá los valores lógicos resultantes de comparar ambos vectores: y %in% x será TRUE, TRUE, FALSE, pues a está en x, d está en x y f, no.

    El conjunto de datos con el que estamos trabajando permite comprobar si los estados de Dakota, Boston y Washington son realmente estados o no. ¿Están los tres elementos de esta lista contenidos en la lista de los 52 estados de EE.UU.? Es decir, c("Boston", "Dakota", "Washington") %in% murders$state. La respuesta es FALSE FALSE TRUE porque ninguno de los dos primeros es un estado, pero Washington sí lo es.

1.8 Manejo básico de datos (basic data wrangling)

Hasta ahora se han manipulado vectores reordenándolos y agrupando sus elementos al indexar el vector de varias maneras (operadores lógicos, datos de otros vectores, etc.). Pero para profundizar en el análisis de datos se necesitan manipular conjuntos de datos enteros, tablas de datos. Para eso es necesario el paquete dplyr.5 Este paquete proporciona las funciones más comunes para el manejo de tablas.

Se pueden añadir columnas (variables) a una tabla o cambiar una columna que ya existe con mutate; se pueden filtrar los datos usando filter o se pueden hacer subconjuntos de datos seleccionando columnas concretas con la función select. Por supuesto, se pueden generar combinaciones de funciones pasando el resultado de una función como parámetro para otra. Para eso se usa la tubería %>%.

1.8.1 Mutate

Esta función modifica columnas en un conjunto de datos o tabla de datos. Necesita dos parámetros: el conjunto de datos con el que va a trabajar y el nombre y valor de la columna que se va a añadir mutate(datos, variable = valor). Tanto la variable como los datos (en el caso de que se obtengan manipulando variables ya existente tal y como pasará a continuación) se extren directamente de la tabla de datos y no del espacio de trabajo, así que es suficiente con el nombre de la variable (population, state, etc.). No es necesario especificar el conjunto de datos una y otra vez (murders$population, murders$state, etc.).

Por lo tanto, para añadir el índice de asesinatos al conjunto de datos como una columna más se usará: murders <- mutate(murders, rate = total/population*100000). Nótese la falta de murders$ delante de cada nombre de variable. Dado que ya se había definido la variable murderrate, también sería válido un código como murders <- mutate(murders, rate = murder_rate).

Ahora, la estructura del conjunto de datos ya incluye esta columna. Se puede comprobar fácilmente con str(murders).

NOTA IMPORTANTE: Esta función ha modificado el objeto murders, que existe en el espacio de trabajo. Si en algún momento se vuelve a cargar este desde el paquete original (library(dslabs)), el objeto se sobrescribirá y se perderán los cambios hechos.

1.8.2 Filter

Por ejemplo, si queremos visualizar solo los datos para los que el índice de asesinatos es menor de 0.71. El primer argumento es la tabla o conjunto de datos y la condición que queremos filtrar, el segundo.

filter(murders, rate <=0.71)

Igual que con mutate, los parámetros son las variables de la tabla directamente, no objetos del espacio de trabajo, por lo que no hay que escribir murders$ constantemente.

1.8.3 Select

La función select permite trabajar con solo el conjunto de variables que se seleccione como segundo argumento (el primero es el nombre de la tabla). En este ejemplo solo hay seis columnas, pero muchas tablas hay cientos. Como ejemplo, crearemos un nuevo objeto con tres columnas (nombre del estado, región a la que pertenece e índice de asesinatos) y luego aplicaremos un filtro:

new_table <- select(murders, state, region, rate)
filter(new_table, rate <= 0.71)

De nuevo, tanto la tabla como sus variables no son objetos del espacio de trabajo.

1.8.4 Las tuberías (pipe): %>%

Tal y como hacen las tuberías, se puede canalizar el resultado de una función para que haga las veces de entrada de datos de la siguiente.

Así, funciones como las anteriores pueden ser escritas en una sola orden. Primero se toma un origen de datos, una tabla murders. Estos datos se «canalizan» hacia una segunda función en la que se seleccionan solo ciertas variables:

murders %>% select(state, region, rate) %>% filter(rate <= 0.71)

El operador de canalización (%>%) entiende que lo que sea que se canaliza es sobre lo que hay que efectuar las operaciones, de ahí que no sea necesario ese primer parámetro que especifica el origen de datos como en los ejemplos de las funciones una por una dados antes. También se evita el objeto intermedio newtable.

Si no se carga la biblioteca dplyr antes, estas funciones (incluida la tubería) no estarán disponibles y se mostrará un error.

1.9 Crear tablas de datos

Para muchos de los análisis que se pueden hacer con el paquete dplyr puede que necesitemos crear una tabla de datos con valores. Una tabla de datos o data frame se puede crear con la función data.frame. En ella se especifican unas variables y valores que asignar a esas variables, como en el ejemplo siguiente:

notas <- data.frame (nombres = c("Juan", "Pedro", "María", "Raquel",
		    "Margarita"),
		    examen_1 = c(2.5, 4, 6.5, 8, 4),
		    examen_2 = c(5.5, 6, 7, 7.5, 3),
			stringsAsFactors = FALSE)
notas

# nombres    examen_1   examen_2
# Juan       2.5        5.5
# Pedro      4          6
# María      6.5        7
# Raquel     8          7.5
# Margarita  4          3

OJO Por defecto, data.frame hace que las entradas alfanuméricas sean tratadas como factores (ver más arriba), por lo que será necesario (si lo es, que puede que en algún caso sea conveniente) añadir el último de los parámetros del código anterior: stringsAsFactors. Si no se incluye, el resultado de class(notas$nombres) será factor. Si se incluye, será character.

1.10 Gráficos básicos

Los gráficos son tremendamente útiles a la hora de inferir algo de los conjuntos de datos de manera intuitiva, mucho más que las listas de números, como pasa con el conjunto de datos del ejemplo, murders.

R es muy potente a la hora de hacer gráficos y mostrar datos.

La función plot hace gráficos de manera muy sencilla.

population_in_millions <- murders$population/10^6
total_gun_murders <- murders$total

plot(population_in_millions, total_gun_murders)

Es decir, plot(datos_en_el_eje_x, datos_en_el_eje_y).

Una vez que se escriba esta orden en una consola de R, R mostrará el gráfico y no parece que haya manera humana de deshacerse de él. En los IDE, como RStudio, hay una sección especialmente dedicada a los gráficos, de manera que es más sencillo gestionarlos (los gráficos). Trabajando en consola, parece que la única manera de cerrar el gráfico es cerrar la ventana que lo contiene.

Creo que hay una manera de guardar el gráfico a disco, pero aún no lo tengo claro.

Los histogramas son muy útiles y se describirán más tarde en el programa, y se consiguen con la función hist, así: hist(murders$rate).

Una manera más detallada de mostrar datos son los diagramas de caja o boxplot ver más información aquí. Un ejemplo de diagrama de caja (boxplot en inglés) es la siguiente: boxplot(rate~regions, data = murders). Esta representación de datos da mucha más información acerca de cómo se distribuyen los datos de una muestra.

2 Programación básica

2.1 Condicionales (lo básico)

2.1.1 if-else

La sentencia condicional más básica es if-else. La fórmula general es la siguiente:

if(condición lógica){
   expresión_si_es_cierta
   } else {
   expresión_alternativa
   }

Un ejemplo con el conjunto de datos de asesinatos por arma de fuego en EE.UU. que hemos usado hasta ahora:

# cargar los datos (por recordar, más que nada)
library(dslabs)
data(murders)

# definir el índice por cada 100,000 habitantes
murder_rate <- murders$total/murders$population*100000

# objeto intermedio para encontrar el menor índice de asesinatos de
# todo el conjunto de datosæ
ind <- which.min(murder_rate)

# sentencia condicional
if(murder_rate[ind]) < 0.5){
   print(murders$state[ind]
   } else {
   print("Ningún estado tiene un índice tan bajo. Vete a Japón,
   pringao")
   }

En esta función condicional se usa un objeto intermedio (ind) que almacena el número de índice (el número de fila) del menor de todos los valores de murderrate. Después, en la condición propiamente dicha, se pregunta si el valor de murderrate que corresponde a esa fila es menor de 0.5. Nótese que ind es un número de fila, mientras que murderrate[ind] es un valor concreto de número de asesinatos por cada 100,000 habitantes, por eso se pregunta por el valor de la variable murderrate filtrada por el número de fila del valor mínimo. Eso devuelve el valor mínimo.

Si la condición es cierta, se muestra en pantalla (print) el valor almacenado en la variable murders$state en la posición (número de fila) correspondiente a ind (el valor mínimo de murderrate).

Si no es cierta (else), se muestra un mensaje. Perfectamente se podría mostrar otra cosa, como el valor máximo, la lista completa o lo que nos dé la gana, esto es solo un ejemplo. No tiene por qué ser un mensaje.

El ejemplo muestra [1] "Vermont", porque el índice de asesinatos en Vermont es murder_rate[ind], o sea, 0.3196. Si en lugar de poner una condición a 0.5, esta se pone a 0.1, el resultado será el mensaje alternativo propuesto: Ningún estado tiene un índice tan bajo, ¡vete a Japón, pringao!

2.1.2 ifelse

La función ifelse está relacionada con la anterior, pero no es lo mismo. Se escribe todo junto, en una sola palabra, y tiene tres parámetros obligatorios: una condición lógica y dos respuestas posibles (si se cumple la condición y si no se cumple).

Un ejemplo que devuelve el recíproco de un número x6 solo si x es positivo es el siguiente:

x <- 0
ifelse(x > 0, 1/x, NA)

La utilidad de esta función reside en el hecho de que funciona con vectores. Efectúa las comprobaciones elemento por elemento, tal y como se ha visto anteriormente en otras funciones. Por ejemplo, si se define a como un vector con números tanto positivos como negativos, la función ifelse comprobará cada uno de los elementos del vector:

a <- c(1, 2, -4, 5)
ifelse(a > 0, 1/a, NA)

#[1] 1, 0.5, NA, 0.2

Se toma un elemento, se comprueba la condición y se devuelve una de las dos respuestas dependiendo de si la condición se cumple o no. Es muy útil (esta función) para sustituir los NA por ceros, por ejemplo. Reemplazar los valores que faltan (NA) por otro valor es un uso muy habitual de esta función. En el paquete dslabs hay otro conjunto de datos de ejemplo con este propósito: naexample.

library(dslabs)
data(na_example)

# La suma de NA se consigue preguntando por los NA contenidos en el
# conjunto de datos y sumándolos.
sum(is.na(na_example))
#[1] 145

# Nuevo vector en el que se sustituyen todos los NA por 0
no_nas <- ifelse(is.na(na_example), 0, na_example)

# Nueva comprobación:
sum(is.na(no_nas))
#[1] 0

En este ejemplo, la función ifelse comprueba si el dato del vector naexample es NA. Si lo es entonces la condición is.na(na_example es TRUE y aplica la primera respuesta, es decir, 0. Si no lo es, aplica la segunda respuesta, es decir, el dato de naexample sin modificar.

2.1.3 any

Esta función comprueba todos los elementos de un vector lógico y devuelve TRUE si cualquiera de los elementos (al menos uno) es TRUE.

z <- c(TRUE, TRUE, FALSE)
y <- c(FALSE, FALSE, FALSE)
any(z)
#[1] TRUE
any(y)
#[1] FALSE

2.1.4 all

Esta función comprueba todos los elementos de un vector lógico y devuelve TRUE si todos los elementos son TRUE.

z <- c(TRUE, TRUE, FALSE)
y <- c(TRUE, TRUE, TRUE)
all(FALSE)
#[1] FALSE
all(y)
#[1] TRUE

2.2 Creación de funciones

A medida que se avanza en el análisis de datos, será necesario repetir la misma operación una y otra vez, como por ejemplo calcular la media de una serie de datos. Se puede hacer a mano, por supuesto, dividiendo la suma de los elementos sum(x) entre la cantidad de elementos length(x). Pero sería mucho mejor escribir una función que lo hiciese, de manera que con escribir media(x) hiciese el cálculo completo.

En realidad no es necesario, porque esa función ya existe, es mean(x), pero otras funciones que aún no existen pueden ser necesarias.

Para ello está la función function, que dice a R que se está definiendo una nueva función, con sus parámetros necesarios u opcionales, configuraciones por defecto, etc.

El uso básico de function es el siguiente:

media <- function(x){
    s <- sum(x)
    n <- length(x)
    s/n
    }

Se almacena en media (ese será el nombre de la nueva función) una nueva función con un parámetro obligatorio entre paréntesis (x) sobre el que se efectúan las operaciones entre llaves: se asigna la suma de los elementos a s, se asigna la longitud del vector a n y se efectúa la división de s entre n. Lo que aparezca en la última línea será la salida de la nueva función, que puede ser un valor a devolver, como en este caso, o una modificación de una tabla, una asignación de valor, el borrado de un dato, etc.

Nota importante: Las asignaciones a variables u objetos que se efectúen dentro de la función no afectan al espacio de trabajo. Es decir, que si hay un objeto s con un valor en el espacio de trabajo porque en algún momento se ha definido s <- 1, este valor no se ve en absoluto afectado por la asignación a s de valores diferentes a 1 (el valor en el ejemplo) durante la operación media. De igual manera, la operación media no tiene en cuenta los valores que s o n pudieran tener en el espacio de trabajo y por tanto no le afectan para nada.

La función general es:

función_nueva <- function(x){
      operaciones a realizar sobre x
      resultado
      }

Puede ser que sean necesarios varios parámetros, en cuyo caso se especifican entre los paréntesis: function(x, y, z){. Por ejemplo, si se piden tres valores de una ecuación cuadrática para dar las soluciones. En ese caso, las operaciones entre llaves se efectuarán sobre los tres elementos aportados por el usuario.

También es posible que los parámetros tengan opciones por defecto. El siguiente ejemplo calcula la media aritmética o geométrica dependiendo de un parámetro a elección del usuario. Este parámetro tiene un valor predefinido que será tomado por defecto:

media <- function(x, aritmetica = TRUE){
    n <- length(x)
    ifelse(aritmetica, sum(x)/n, prod(x)^(1/n))
    }

En la última línea, una sentencia ifelse hace que, dependiendo de si el parámetro aritmetica existe o no (es decir, si es 1 o 0, si es TRUE o FALSE) devuelve la media aritmética (es la primera opción de ifelse, la que se efectúa si la condición es cierta), y este es el comportamiento por defecto, pues por defecto el parámetro aritmetica se ha establecido como TRUE. Si, por el contrario, hay una definición explícita del aritmetica como FALSE, entonces ifelse escoge la segunda respuesta y devuelve la media geométrica.

2.3 Bucles (bucles for)

Los bucles se crean para repetir la misma operación una y otra vez sobre un conjunto de datos en el que se varía el dato sobre el que se opera en cada una de las iteraciones del bucle.

Por ejemplo, para calcular la suma de una progresión. La progresión de los números naturales desde 1 hasta n 1+2+3+4+5+…+/n/ se puede calcular mediante la fórmula [/n/(n/+1)1]/2. Para comprobar si la fórmula es correcta, se pueden efectual las sumas desde 1 hasta /n solo escogiendo un número n.

Una función que haga el cálculo es sencillo después de ver el apartado anterior sobre creación de funciones:

suma_n <- function(n){
     x <- 1:n
     sum(x)
     }

La función pide un parámetro y asigna al objeto x un vector formado por los números naturales desde 1 hasta n. A continuación calcula la suma de los elementos de ese vector.

Ahora bien, si lo que necesitamos no es conocer la suma hasta un n concreto sino cómo evoluciona esta suma desde n=1 hasta n=25 es necesario hacer el mismo cálculo 25 veces, cada vez con un valor diferente de n. Para eso vale un bucle for.

La forma general es:7

for(i in rango_de_valores){
    operaciones a realizar sobre i a medida que varía el valor de i
    }

Un bucle extremadamente sencillo (e inutil) es el siguiente:

for(i in 1:5){
   print(i)
   }

Da como resultado los números de uno a cinco. Lo que hace es asignar valores a i dentro de los definidos entre paréntesis (de uno a cinco) y luego ejecuta el código sobre ese valor de i. En este caso es solamente mostrar i.

Nota importante: el último valor de i (o cualquier otra variable que se use en el bucle for) se mantiene después del cálculo. Con el ejemplo anterior, el código i devolverá 5, porque es el último valor de i con el que se ejecutó el bucle.

Para el ejemplo propuesto, la suma de los primeros n números naturales, se necesita primero crear un vector que almacene los n resultados de las n sumas desde 1 hasta n: m <- 25 para definir las sumas desde 1 hasta 25, y s_n <- vector(length = m) para crear un vector vacío de esa longitud.

Ahora viene el bucle en sí. Nótese que suma_n(n) es la función creada anteriormente.

for(n in 1:m){
   s_n[n] <- suma_n(n)
   }

Este código define un valor cambiante de n desde 1 hasta m y asigna el valor de la función suman sobre ese valor de n. La función suman definida anteriormente calcula la suma desde 1 hasta n y la almacena en la posición n en el vector sn. Ese vector tiene una longitud m para poder almacenar tantos valores de n como se definan en la variable m.

Un gráfico mostrará hasta qué punto la secuencia lograda es correcta: plot(1:m, s_n). Además, se puede superponer la funcíon lines para mostrar el resultado de la fórmula y comparar la superposición de ambos resultados tal que así: lines(n, n*(n+1)/2). Sin embargo, en consola al menos no funciona (plot muestra el gráfico, pero lines da una pantalla en blanco. Hay que probar en RStudio u otro IDE.

2.4 Otras funciones útiles

Este curso no las cubrirá, pero funciones útiles que es provechoso aprender a utilizar son:

  • Funciones de la familia apply:
    • apply
    • sapply
    • tapply
    • mapply
  • split
  • cut
  • quantile
  • reduce
  • identical
  • unique

Lo que se aprenda de ellas irá apareciendo aquí.

3 Visualización de datos

3.1 Mostrar y guardar los gráficos

Trabajando en la consola de R los gráficos salen directamente a la pantalla y no se puede hacer nada con ellos, salvo matar la ventana y cerrarlos. Para guardar las imágenes hay que hacer lo siguiente:

  1. Llamar a la función correspondiente a la imagen que se quiere crear.
  2. Ejecutar el gráfico, que no aparece.
  3. Liberar el control y devolverlo a la pantalla.

Hay varias funciones para cada uno de los tipos de gráfico: jpeg(), png(), bmp() y tiff() para mapas de bits.

jpeg(file = "gráfico_1.jpeg")
plot(murders$region)
dev.off()

Se puede especificar la ruta completa para no guardar la imagen en el directorio actual. También se puede especificar el tamaño de la imagen con width y height. Por defecto es 480x480. Por ejemplo, guardando como png.

png(file = "gráfico_2.png",
    width = 640,
    height = 640
    )
plot(murders$region)
    dev.off()

Argumentos posibles son las unidades en las que se mide el tamaño (mm, cm, in) y la resolución (en ppp), como en

png(file = "gráfico_3.png",
    width = 45
    height = 30
    units = "mm"
    res = 300
    )

Otra opción es guardar como pdf con pdf(file = "archivo.pdf") o como PostScript con postscript(file = "archivo.ps").

3.2 Introducción a la visualización de datos y distribuciones

3.2.1 Introducción a la visualización de datos

Corresponde al capítulo 6 del libro

Normalmente los datos en bruto no son muy útiles a la hora de derivar ideas de una tabla o un vector. Los números por sí solos no dicen mucho a la mayoría de las personas. Sin embargo, un gráfico puede ofrecer mucha información fácilmente comprensible de un vistazo, de manera que se pueden extraer conclusiones de manera más sencilla que mirando tablas. En ocasiones, una imagen es tan demostrativa que no hace falta ningún análisis más.

Últimamente, cada vez hay más datos disponibles y herramientas informáticas para analizarlos. Incluso en el mundo del periodismo se aprecia un aumento de los gráficos basados en datos para apoyar afirmaciones.

La visualización correcta de los datos es la herramienta más poderosa del análisis exploratorio. Y este es a su vez la parte más importante del análisis de datos.

3.2.2 Introducción a la distribución normal

Corresponde al capítulo 8 del libro. Es habitual resumir conjuntos de datos numéricos con un solo valor (la media) y, en ocasiones, con dos valores (la media y la desviación típica): el rendimiento de una clase con el valor medio de la nota de un examen (5.2 o 5.2±0.4).

En ocasiones estos dos números son suficientes para entender los datos. Las técnicas de visualización de datos ayudan a ver cuándo estos resúmenes son adecuados y proporcionan herramientas y alternativas cuando no lo son.

Lo más básico para analizar un conjunto numérico es la distribución. Una vez que un vector ha sido resumido como una distribución, hay técnicas que permiten extrer información muy útil de ese vector.

  1. Tipos de datos

    En un conjunto de datos (data set), hay datos catalogados como registros (cada una de las mediciones hechas, como cada persona que responde una encuesta, las filas de la tabla) y como variables (cada uno de los datos que cada registro genera, como las preguntas de una encuesta, las columnas de la tabla).

    Las variables pueden ser de dos tipos:

    • categóricas
    • numéricas

    Cuando una variable tiene un número muy limitado de valores, se dice que es categórica. Ejemplo: sexo (hombre, mujer), regiones en EE.UU. (South, West, NorthEast, Central).

    Algunos datos categóricos pueden ser ordenados aunque no sean números, como los grados de temperatura de un producto: frío, templado, caliente. Pueden ser conocidos como datos ordinales (se pueden ordenar).

    Las variables numéricas pueden ser continuas o discretas. Una variable continua es una que puede tomar cualquier valor, como la altura de una persona (1.82, 1.772). Los valores que puede tomar dependen únicamente de la precisión con que se mida. Podría ser considerada categórica solo si se toman los valores alto, mediano y bajo.

    Una variable discreta es una que solo puede tomar un cierto número de valores a números redondos, como la altura en centímetros (182, 177) o una población (233,845 personas). Según algunos textos, se pueden considerar ordinales también, pero puede ser confuso. Es mejor ceñirse a la nomenclatura de variables ordinales como aquellas que tienen un rango muy pequeño de valores. Por ejemplo, el número de paquetes de cigarrillos que una persona fuma al día puede ser considerada ordinal (valores 0, 1, 2 o 3), pero el número de cigarrillos será una variable numérica (de 0 a 60).

    Se puede conocer el nombre de las variables de un conjunto de datos con la función names.

    install.package("dslabs")
    library(dslabs)
    data(heights)
    names(heights)
    
    # [1] "sex"  "height"
    

    Se pueden conocer los valores que puede tomar una variable observando los primeros registros de una tabla con la función head:

    library(dslabs)
    data(heights)
    head(heights)
    
    #   sex  height
    #  1   Male     75
    #  2   Male     70
    #  3   Male     68
    #  4   Male     74
    #  5   Male     61
    #  6 Female     65
    

    COMPLETAR

  2. Un primer ejemplo

    Un problema artificial pero ilustrativo para entender los conceptos detrás de la idea de distribución: describir la estatura de un grupo de personas a ET (un extraterrestre que nunca ha visto un humano). Lo primero es recolectar los datos. Eso lo consideramos hecho con el conjunto de datos heights ya incluido en el paquete dslabs.

    El modo más sencillo es enviar el conjunto completo de datos con sus 1,050 entradas. Pero hay maneras más efectivas.

    La distribución es el resumen estadístico básico de una lista de objetos o números, es una descripción compacta de la lista. Esto es muy sencillo con datos en categorías, como la lista de sexos en heights: male y female. Así que el resumen podría ser simplemente la proporción de ambos. El código prop.table(table(heights$sex)) dará las proporciones en tanto por uno, o sea Female será 0.227 y Male será 0.773.8 Y no hay más.

    Para categorías más complejas, un gráfico de barras describe las proporciones de cada categoría. Son ejemplos de cómo convertir un vector en un elemento visual que resuma sus características.

    Con datos numéricos es algo más complejo, pues como cada elemento es único (o casi, la frecuencia de repetición suele ser muy baja), mostrar la frecuencia de aparición no es muy efectivo. Sería mejor definir una función que mostrase la proporción de datos por debajo de un valor A para todos los posibles valores de A (esto es una distribución acumulativa o CDF por sus siglas en inglés (Cumulative Distribution Function). Es decir, cuántos valores hay de A entre el inicio y A. Comprobando el valor de la curva en x = 72 se ve que el valor es 0.839, es decir, el 83.9% de los elementos de la curva tiene una altura igual o menor a 72 pulgadas. ¿Cuántos alumnos tienen una altura entre 66 y 72 pulgadas? Son el 83.9% menos el 13.8% = 70.1%. Es por tanto muy útil para decidir la probabilidad de que una medida tenga un valor determinado (un 70% de que el alumno tenga una altura entre 66 y 72 pulgadas)

    ecdf-1.png

    Figure 1: Gráfico de distribución acumulativa

    Sin embargo, no es un gráfico muy popular porque obvia datos interesantes como el valor en el que se centra la lista de datos, si es simétrica o no respecto a ese centro, etc. Los datos están, pero no son fácilmente identificables.

    Este gráfico se puede calcular con el siguiente código:

    # Dado un conjunto de datos llamado *datos*
    # se define el rango de valores que abarcan los datos
    a <- seq(min(datos), max(datos), length = 100)
    
    # se calcula la probabilidad para un valor
    cdf_function <- function(x){
        mean(datos <= x)
    }
    cdf_values <- sapply(a, cdf_function)
    plot(a, cdf_values)
    
  3. Histogramas

    Una manera correcta y muy exacta de representar datos es con una gráfico de distribución acumulada como el anterior, pero no es muy práctico en la vida real. Lo más esclarecedor es un histograma. En un histograma se divide una secuencia continua de datos (como las alturas individuales de un grupo de personas) en tramos discretos en los que se cuenta cada una de las apariciones de datos (por ejemplo, entre 176 y 178 centímetros hay 48 personas, entre 178 y 180 hay 30, etc.).

    histo.png

    Figure 2: Histograma

    Es más sencillo ver cómo se distribuyen los datos en la muestra usando un histograma: rango de datos, rango en el que se agrupa la mayoría de elementos, simetría de los datos, etc.

    Pero se pierde algo de información porque dentro de cada intervalo considerado (el ancho de las barras) todas las medidas son consideradas iguales. En el gráfico del ejemplo más arriba cada barra representa intervalos de 5°C. Todas las medidas entre 70 y 75, por ejemplo, son consideradas iguales aunque no lo sean en realidad. Las implicaciones prácticas de esta normalización pueden ser irrelevantes dependiendo del caso.

  4. Curva de densidad

    Una curva de densidad (Smooth Density Plot) es mejor para una interpretación rápida de la distribución de datos. Es una estimación de cómo se vería un histograma con una cantidad extremadamente elevada de observaciones y, por tanto, de barras. Es como aumentar enormemente la resolución con la que se puede ver un histograma.

    smoot-densi.png

    Figure 3: Curva de densidad

    En este tipo de gráficos, el eje Y muestra la densidad de aparición de un valor determinado en el conjunto de datos. Es decir, para conocer la frecuencia de aparición de un cierto valor dentro de un intervalo, hay que conocer el porcentaje del área bajo la curva (la integral de la función) dentro del dicho intervalo en el eje X.

    Es útil para comparar dos muestras. Pero tiene como contrapartida que es una estimación que puede ser ajustada y por tanto se presta a interpretación. Su utilidad depende en parte de la habilidad del analista.

    Su origen es la asunción de que la lista de datos (en el ejemplo es la lista de alturas de alumnos varones en pulgadas contenida en el paquete dslabs) proviene de una lista hipotética que contiene todas las alturas de todos los estudiantes varones de todo el mundo medidas con asombrosa precisión. Y esta es la lista que queremos mostrar a ET, que proporcionará valores más precisos y generales que la lista que tenemos. Pero para eso hay que asumir algo. La precisión con la que se miden los datos y el número enorme de datos permite hacer un histograma con las barras extremadamente delgadas, de manera que no hay saltos entre barras (de ahí lo de smooth 'suave').

    Pero no tenemos infinitas ni un millón de observaciones sino 812 en nuestra lista. Por eso se hace un histograma contando frecuencias de valores en lugar de contar valores. Se dibuja entonces una curva que conecte los picos superiores de las barras del histograma. Cómo se ajusta la curva a la forma del histograma original es algo, como se dijo, que puede ser ajustado con un parámetro en la función ggplot, y depende del analista. En nuestro caso, se puede asumir que las proporciones de medidas contiguas serán similares, de manera que el gráfico tenderá a ser más suave que en otros casos.

    Para entender una curva de distribución como es debido hay que entender cómo se muestran los valores en el eje Y. Como se dijo, no es un gráfico basado en un conteo de elementos sino en la frecuencia de aparición. Por lo tanto, la curva entera tiene una frecuencia de aparición (o probabilidad de aparición de una medida determinada) de 1, es decir, el área bajo la curva es 1. Pero, ¿qué quieren decir los valores de Y?. Son, efectivamente, la frecuencia de aparición de cada medida, y dan esta información si y solo si el intervalo medido en el eje X es 1. Si no lo es, lo que muestra la frecuencia de aparición de este intervalo de valores es el área bajo esa porción de la curva.

    Dicho de otra manera, la proporción de elementos en un determinado rango de medidas en el eje X se calcula con la proporción del área de la curva delimitada por los valores en el intervalo que se quiere medir. Si es un 33% del área total bajo la curva, entonces un 33% de los elementos están contenidos en ese intervalo.

    Es muy útil para comparar dos muestras diferentes, por ejemplo los conjuntos de datos Male y Female en el ejemplo de las alturas mostradas a ET. El paquete ggplot tiene opciones muy potentes para dibujar curvas con áreas coloreadas de manera diferente. En la imagen, estos son los conjuntos que se comparan:

    smoot-compa.png

    Figure 4: Comparación de muestras

  5. La distribución normal

    Tanto los histogramas como la curva de densidad son excelentes resúmenes de una distribución de datos. Para resumir y explicar los datos aún mejor y de manera más resumida se usan la media y la desviación típica (así que se puede resumir la distribución con solo dos números). Esto es gracias a la distribución normal.

    La distribución normal también se conoce como «Campana de Gauss» o «Curva de Bell», y es uno de los conceptos matemáticos más usados de la historia porque es una distribución que se repite en multitud de ocasiones de manera más o menos exacta (en resultados de apuestas, estatura de un grupo, presion arterial, medición de errores experimentales, etc.).

    La distribución normal no se basa en datos empíricos sino en una fórmula matemática. Esta formula se controla en el intervalo entre a y b con tan solo dos parámetros, m y s (la media y la desviación típica). El resto son constantes conocidas (pi y e).

    Lo más importante es que la distribución es simétrica alrededor de la media (igual hacia un lado que hacia el otro) y el 95% de los valores están dentro de ±2 desviaciones típicas de la media. Este hecho, que la formula y por tanto la curva puede ser descrita con solo dos parámetros (m y s) hace que cuando un conjunto de datos es aproximadamente como la distribución normal, toda la información necesaria para describir este conjunto de datos y las inferencias que se pueden extraer se pueden codificar con esos dos parámetros, la media y la desviación típica que (ahora sí) se extraen empíricamente de la muestra.

    La media es la suma de los elementos dividida por el número de elementos. En un vector x la media es media <- sum(x) / length(x), aunque existe una función ya definida para ello: mean(x).

    La desviación típica es la raíz cuadrada de la suma de los cuadrados de las diferencias entre los valores y la media, todo dividido por el número de elementos. Es decir, es la media (cuadrado y raíz cuadrada se anulan pero se usan para eliminar los negativos) de la distancia entre los valores y la media. En R se calcula con:

    desv_tip <- sqrt(sum((x-media)^2) / length(x))`
    

    Es importante saberlo porque la función definida en R no divide por length(x) sino por length(x) - 1.9 Cuando el número de elementos es grande no tiene ninguna importancia un elemento arriba o abajo, así que se puede usar la función sd sin problema. Cuando son pocos, esta diferencia influirá en el resultado.

    Calculemos la media y la desviación típica para el conjunto de datos heights, pero solo los varones. Se guarda en x el vector de estaturas para los varones:

    library(dslabs)
    data(heights)
    index <- heights$sex == "Male"
    x <- heights$height[index]
    # se da una vuelta para mostrar los resultados como una sencilla
    # tabla, pero valdría con las dos funciones tal cual.
    media <- mean(x)
    desv_tip <- sd(x)
    c(media=media, desv_tip=desv_tip)
    
    #   media   desv_tip
    #   69.31       3.61
    

    Un gráfico de densidad muestra que la forma de la curva es prácticamente igual a la distribución normal calculada con la media de 69.31 y la desviación típica de 3.61.

    1. Unidades estándar (estándar units)

      Es conveniente hablar de unidades estándar cuando se trata de conjuntos de datos que se aproximan bastante a la distribución normal. Una unidad estándar informa de cuántas desviaciones típicas se separa un dato contando desde la media. Se calcula como

      z = (x - m)/s

      De esta manera, z se puede usar como medida de cómo de apartado está un valor de la media y, por tanto, cuán extraño es en la muestra. Si z es mayor de 2 o menor de -2, se puede considerar que la estatura es alta o baja (no está dentro de los valores que tiene el 95% de la muestra, que son ±2 desviaciones típicas desde la media -ver arriba-). Si los valores de z son más extremos, pues más alejado el elemento de la media. Cuando el valor de z se acerca a 0, el valor está dentro de la media.

      Otra ventaja es que estos valores de ±2 s alrededor de la media no tienen nada que ver con las unidades originales (pulgadas, mmHg, °C, kg, etc.)

      En R se calcula z con la función scale, así: z <- scale(x). Se puede calcular cuántos alumnos de la muestra están dentro de ±2 desviaciones típicas calculando el valor absoluto de z < 2 y hallando su media (luego lo explico): mean(abs(z) < 2). El valor resultante es 0.9495, ~0.95, o sea, el 95% de las estaturas están dentro del rango de ±2 desviaciones típicas.

      Explicación: en este código se usa la media de los valores absolutos porque el resultado de lo que está entre paréntesis es un vector lógico en el que hay valores TRUE cuando el valor absoluto de z es menor de 2 y valores FALSE cuando es 2 o mayor. Dado que ambos se asimilan a 1 y 0, respectivamente, la media muestra la proporción de unos respecto al total de elementos (suma de unos dividida por el número total de elementos). Si todos fuesen TRUE, la media sería 1.

      De manera que se pueden predecir las proporciones de elementos dentro de un rango sin ni siquiera mirar los datos.

      Sin embargo, para confirmar la exactitud de las observaciones se deben explorar otros indicadores, como los cuantiles. Pero antes, algunos datos más sobre la distribución normal.

  6. La función acumulativa normal y pnorm

    Ya se ha hablado de la distribución normal como aproximación súper útil para muchas distribuciones del mundo real.

    La función aculumativa de la distribución normal viene dada por una fórmula matemática, igual que la curva de Bell, que en R se obtiene con la función pnorm(). Decimos que un valor aleatorio a está normalmente distribuido con media m y desviación típica s si su distribución de probabilidad (ver la función acumulativa o eCDF más arriba) está definida por F(a) = pnorm(a, m, s).

    Otra vez, esto es sumamente útil porque para calcular la probabilidad de que una observación esté contenida en un rango no hace falta todo el conjunto de datos, solo la media y la desviación típica.

    Por ejemplo, para saber la probabilidad de que un estudiante sea más alto que 70.5 pulgadas habrá que calcular 1 - pnorm(70.5, mean(x), sd(x)), donde x es el vector obtenido anteriormente de las alturas de solo los varones. Sin restar de 1 el resultado será la probabilidad de que la altura sea menor de 70.5.

    La distribución normal se calcula matemáticamente solo a partir de la media y la desviación típica, y está definida solo para variables continuas, no tiene sentido en variables discretas. Esto es importante, porque no se puede «preguntar» a los datos acerca de la probabilidad de que la medida sea 70 pulgadas; la probabilidad de que sea 70 pulgadas exactamente es la misma que cualquier otra medida. Hay varias medidas de 70 pulgadas y solo una de 69.6850393700787, pero esto es solo porque la mayoría de los encuestados para este conjunto de datos conoce su altura en pulgadas (porque son estadounidenses, creo) y redondea al valor entero más cercano (seguro que redondean hacia arriba, nadie quiere ser bajito). Pero este valor tan exacto es 177 cm convertido a pulgadas (177/2.54). De manera que, medidos con infinita precisión, todos y cada uno de los valores de la serie de datos serían únicos.

    De lo anterior se deduce que lo que sí tiene sentido es preguntar, en cambio, cuál es la probabilidad de que la medida esté entre 69 y 70 pulgadas. Es decir, las probabilidades se definen para intervalos.

    En casos como heights, donde los datos están redondeados, la aproximación normal es particularmente útil cuando el intervalo incluye uno de los valores después del redondeo, por ejemplo, el intervalo 69.5 a 70.5

    # probabilidades en los datos reales con rangos de 1 pulgada que
    # contienen un valor entero:
    mean(x <= 68.5) - mean(x <= 67.5)
    mean(x <= 69.5) - mean(x <= 68.5)
    mean(x <= 70.5) - mean(x <= 69.5)
    
    # probabilidades calculadas en la distribución normal concuerdan:
    pnorm(68.5, mean(x), sd(x)) - pnorm(67.5, mean(x), sd(x))
    pnorm(69.5, mean(x), sd(x)) - pnorm(68.5, mean(x), sd(x))
    pnorm(70.5, mean(x), sd(x)) - pnorm(69.5, mean(x), sd(x))
    

    Ejecutando el código anterior se ve que las probabilidades en el conjunto de datos reales no están muy alejadas del conjunto normal definido con las mismas media y desviación típica.

    Sin embargo, el siguiente código hace lo mismo para un intervalo que no incluye un valor entero (70.1 y 70.9). La cosa cambia. El valor devuelto (lo probabilidad de encontrar una medida en el intervalo) es 0.0221 para el conjunto de datos, mientras que la distribución normal espera una probabilidad de 0.0836:

    mean(x <= 70.9) - mean(x <= 70.1)
    pnorm(70.9, mean(x), sd(x)) - pnorm(70.1, mean(x), sd(x))
    

    Esto suele llamarse discretización. Aunque la distribución es continua, los valores tienden a agruparse en valores discretos debido al redondeo. Aún así, la distribución normal como aproximación es muy útil, siempre y cuando se tenga esto en cuenta.

    Ejemplos de cómo usar pnorm en un conjunto de datos, en este caso para estimar la proporción de datos entre 69 y 72 pulgadas en el conjunto heights:

    library(dslabs)
    data(heights)
    x <- heights$height[heights$sex=="Male"]
    # Calcular la media (avg) y la desviación típica (stdev) del conjunto
    # de datos
    avg <- mean(x)
    stdev <- sd(x)
    # Calcular la proporción sin el conjunto de datos, solo con los
    # valores de la media y la desv.típ. que los definen:
    pnorm(72, avg, stdev) - pnorm(69, avg, stdev)
    

    Sin embargo, cuando nos acercamos a los bordes de la distribución, puede haber valores en el conjunto de datos donde la distribución normal no los espera, de manera que las diferencias empiezan a ser grandes: el siguiente código muestra que cuando se calcula la proporción de valores entre 79 y 81 pulgadas, el valor real es 1.6 veces el valor normal:

    library(dslabs)
    data(heights)
    x <- heights$height[heights$sex == "Male"]
    exact <- mean(x > 79 & x <= 81)
    approx <- pnorm(81, mean(x), sd(x)) - pnorm(79, mean(x), sd(x))
    exact / approx
    

    Como se dijo antes, pnorm calcula la proporción de valores menores o iguales a un valor dado x. Para saber la proporción de valores mayores hay que restar de 1: 1 - pnorm(x, mu, sigma).

    Cuando se usan proporciones para calcular poblaciones (como en ¿Cuántos hombres en el mundo son más altos que x si la proporción que muestra la función pnorm es 0.04?). En este ejemplo, hay que multiplicar la proporción por el número de hombres, pero no se pueden medir decimales, así que hay que redondear con round:

    # Proporción *p* para el valor 84 en una distribución normal de media *mu* 69 y desviación típica *sigma* 3:
    p <- 1 - pnorm(84, 69, 3)
    
    # Redondeo para aplicar la proporción a una población de 10^9
    # elementos (población de varones entre 18 y 40 años):
    round(p * 10^9)
    

    Eso son 287 hombres (entre 18 y 40 años) más altos de 84 pulgadas (7 pies o, lo que es lo mismo, 213 cm.). Curiosidad: Si hay 10 jugadores de la NBA más altos de 213 cm, ¿qué proporción de la población mundial tan alta está en la NBA? Pues 10 / round(p * 10^9), es decir, un 3.48%.

    Aparte de esta curiosidad, ahora se repite el cálculo para la altura de Lebron James, 6 pies y 8 pulgadas o 203 cm, sabiendo que hay 150 jugadores en lugar de 10 que son tan altos como él. Hay un 0.12% de hombres de esa altura en el mundo jugando en la NBA. Conclusión: es más fácil entrar en la NBA si eres muy alto. FALSO. Como se vio antes, la distribución normal tiende a subestimar las medidas en los extremos, o sea que es muy probable que haya más hombres de 7 pies o más en el mundo de los que espera la distribución normal.

3.2.3 Cuantiles, percentiles y diagramas de cajas (y bigotes)

Los cuantiles son puntos de corte que dividen un conjunto de datos en intervalos de una probabilidad dada. El cuantil q-ésimo es un valor del que el q% de las observaciones es igual o menor. El cuantil q se puede calcular con la función quantile(data, q). El valor q debe proporcionarse en tanto por uno, es decir, el 78% de las observaciones es igual o menor a 71.11 pulgadas, de manera que quantile(heights$height, 0.78) devolverá 71.11 como resultado. Valores mayores que 1 darán error.

Por lo tanto, el cuantil 0 es el valor mínimo, el cuantil 0.5 es la mediana y el cuantil 1 es el valor máximo (el 100% de las observaciones son iguales o menores a ese valor).

Dentro de los cuantiles (que es la fórmula general) se pueden encontrar varias maneras de agrupar los datos. Las más usadas son:

  • Los cuartiles, que dividen a la distribución en cuatro partes (corresponden a los cuantiles 0,25; 0,50 y 0,75);
  • Los quintiles, que dividen a la distribución en cinco partes (corresponden a los cuantiles 0,20; 0,40; 0,60 y 0,80);
  • Los deciles, que dividen a la distribución en diez partes;
  • Los percentiles, que dividen a la distribución en cien partes.

Los percentiles son los cuantiles que dividen un conjunto de datos en 100 intervalos, cada uno con un 1% de probabilidad. Si se define una secuencia desde 0.01 (= 1%) hasta 0.99 en intervalos de 0.01 y se calculan los cuantiles correspondientes a cada elemento de esa serie se tendrá una lista completa de los percentiles de un conjunto de datos. Ejemplo:

p <- seq(0.01, 0.99, 0.01)
quantile(data, p)

Los cuartiles son lo mismo pero para cuantiles 0.25, 0.50 y 0.75. También conocidos como primer cuartil, mediana y tercer cuartil. Una visión general de los cuartiles se puede lograr con la función summary: summary(heights$height) devolverá el valor mínimo (el cuantil 0), el primer cuartil, la mediana, la media, el tercer cuartil y el valor máximo.

  1. Encontrar cuantiles con qnorm

    La función qnorm da el valor teórico del cuantil de probabilidad p de una distribución de media mu y desviación típica sigma. Por defecto, mu es 0 y sigma es 1, por lo que si se escribe sin argumentos qnorm() da el valor del cuantil de una distribución normal teórica estándar. En efecto, qnorm(0.5) devuelve 0, porque en una distribución normal estándar la media mu coincide con la mediana (cuantil 0.5). Es decir, hay tantos valores por encima de la media como por debajo. Para especificar la media y la desviación típica se hace por este orden: qnorm(p, mu, sigma).

    Si pnorm ofrece la probabilidad de que un valor sea menor o igual a z, qnorm es la función inversa y ofrece el valor que tiene una probabilidad p.

    Como qnorm devuelve los valores de los cuantiles teóricos, se puede usar para cerciorarse de si la distribución de valores del conjunto de datos es aproximadamente normal o no. Un ejemplo: conociendo mu y sigma para un conjunto de datos, como el de las alturas de alumnos, los cuantiles teóricos se pueden calcular de la siguiente manera:

    #Elaborar una secuencia desde 1% a 99% en tanto por uno, de uno en uno y asignarla al vector *p*
    p <- seq(0.01, 0.99, 0.01)
    
    #Calcular los cuantiles almacenados en *p* de una distribución
    #(teórica) con *mu* y *sigma* como valores definitorios.
    cuantiles_teoricos <- qnorm(p, 63, 3)
    
  2. Diagramas q-q o cuantil-cuantil y diagramas de cajas

    Para saber si la distribución normal es una buena aproximación a la distribución de los datos en el conjunto estudiado se puede hacer una gráfica comparando los cuantiles teóricos con los observados. Si ambos coinciden (más o menos), se puede usar una distribución normal de media mu y desviación típica sigma como resumen del conjunto de datos completo. Insisto: este es un método para hacer la comprobación de la validez de la distribución normal como resumen del conjunto de datos; si es válida, el conjunto de datos se resumirá con solo dos valores (la media y la desviación típica).

    Antes de crear un gráfico, se pueden comparar dos conjuntos de datos con una sencilla tabla en la que poner lado a lado varios percentiles de ambos conjuntos de datos. Por ejemplo, comparar los percentiles 10º, 30º, 50º, 70º y 90º de las alturas de alumnos varones y mujeres.

    # Cargar los datos
    library(dslabs)
    data(heights)
    
    # Asociar los datos filtrados a sendos vectores
    male <- heights$height[heights$sex == "Male"]
    female <- heights$height[heights$sex == "Female"]
    
    # Definir un vector con los percentiles que se van a comparar
    percntls <- seq(0.1, 0.9, 0.2)  # de 10% a 90% en saltos de 20
    
    # Calcular los percentiles y asociarlos a vectores
    female_percentiles <- quantile(female, percntls)
    male_percentiles <- quantile(male, percntls)
    
    # Crear un objeto que contenga un conjunto de datos (data frame)
    df <- data.frame(female = female_percentiles, male = male_percentiles)
    
    # Comprobar el resultado
    df
    

    Ahora sí, el código necesario para crear un gráfico que permita la comparación es el siguiente:

    # Cargar los datos
    library(dslabs)
    data(heights)
    
    # Asociar datos a un vector y filtrar
    index <- heights$sex == "Male"
    x <- heights$height[index]
    z <- scale(x)
    
    # calcular los cuantiles teóricos y reales
    p <- seq(0.05, 0.95, 0.05)
    observed_quantiles <- quantile(x, p)
    theoretical_quantiles <- qnorm(p, mean = mean(x), sd = sd(x))
    
    # Dibujar el gráfico
    plot(theoretical_quantiles, observed_quantiles)
    abline(0, 1)
    

    El código es más sencillo si se usan unidades estándar (número de desviaciones típicas desde la media), ya que no hay que calcular la media y desviación típica para el valor teórico de referencia:

    # Gráfico con valores escalados (en desviaciones típicas desde la
    # media)
    observed_quantiles <- quantile(z, p)
    theoretical_quantiles <- qnorm(p)
    plot(theoretical_quantiles, observed_quantiles)
    abline(0, 1)
    

    Si el gráfico se dibuja en una consola de R (no en un IDE), se puede volver a la ventana de terminal e introducir la orden abline(0, 1) mientras el gráfico está en otra ventana. La línea será dibujada en el gráfico directamente.

3.2.4 Análisis exploratorio

Haciendo lo mismo que se hizo con las alturas de los estudiantes varones, pero ahora con las alturas de las mujeres, se ve que hay variaciones importantes respecto a la distribución normal: hay un aplanamiento de la curva en torno a un valor superior a la media y en un gráfico q-q se ve que hay más mujeres muy altas de lo esperado.

Esto nos puede ayudar a detectar problemas en el conjunto de datos. Por ejemplo, en el formulario on-line en el que los estudiantes introducían su altura, la opción por defecto era Sexo: mujer y hubo chicos que detallaron su altura pero no cambiaron la opción en cuanto al sexo.

Los cálculos básicos para saber algo de los datos disponibles son: 1. La media mean(x). 2. La mediana median(x). 3. La desviación típica sd(x). 4. La desviación absoluta respecto a la mediana10 mad(x).

Tanto la media como la desviación típica son muy sensibles a las variaciones de valores fuera del rango típico (outliers), por ejemplo errores de tecleo, como comerse puntos decimales. La mediana y la desviación absoluta respecto a la mediana son mucho menos sensibles y se denominan medidas robustas.

Estos valores que se salen del rango claramente son detectables con un diagrama de cajas, un histograma o un gráfico q-q.

Se puede escribir una función para comprobar cómo estos valores fuera de rango o outliers afectan a la media de una serie de datos. Se toma una serie de datos, se varía uno de ellos y luego se comprueba la media, por ejemplo con un valor verídico pero sin coma decimal (es decir, multiplicado por 10).

# Cargar los datos
library(dslabs)
data(heights)
index <- heights$sex == "Male"
x <- heights$height[index]

# Escribir la función
error_avg <- function(k){
  index <- heights$sex == "Male"
  x_err <- heights$height[index]
  x_err[1] <- k
  mean(x_err)

# Comprobar el resultado (69.5 * 10)
error_avg(695) - mean(x)

En este caso, la media sube más de una pulgada y media con solo un valor erróneo en todo el conjunto.

3.3 Introducción a ggplot2

ggplot2 es relativamente sencillo de usar para principiantes, comparado con otras herramientas para crear gráficos. Se instala con el conjunto de paquetes tidyverse: library(tidyverse) o puede cargarse ggplot2 en exclusiva con library(ggplot2). Antes de cargar las bibliotecas puede ser necesario instalarlas: install.packages("tidyverse"). Al hacerlo, atención a los mensajes de R: en mi caso he necesitado instalar primero tidyr y broom.

ggplot está diseñado para trabajar exclusivamente con conjuntos de datos ordenados en tablas, donde las filas son observaciones y las columnas son las distintas variables.

3.3.1 Componentes del gráfico

Hay varias partes en las que se puede dividir un gráfico para ver más fácilmente cómo se construye y qué órdenes se incluyen en la función ggplot2 para lograr un resultado conveniente.

  1. Datos: el conjunto de datos incluido en el gráfico.
  2. Geometría: diagrama de dispersión, diagrama de barras, histograma, etc.
  3. Aspecto: valores de los ejes, colores de puntos, barras o líneas, etc.
  4. Otros componentes como escala de los ejes, estilo, leyenda, títulos, etc.

3.3.2 Crear un nuevo gráfico

Lo primero que hay que hacer es definir un «objeto ggplot» con la función ggplot. Esta función asocia un conjunto de datos con el gráfico que se va a crear (el componente número 1, datos, visto en el apartado anterior). Pero aún no hay nada dibujado, tan solo un objeto definido y los datos correspondientes asociados a él. En realidad, el código sí dibuja un gráfico, pero lo único que se ve es un fondo gris, falta la geometría, el aspecto y otros componentes necesarios.

Para definir el objeto, las líneas de código siguientes son equivalentes:

ggplot(data = murders)
murders %>% ggplot()

p <- ggplot(data = murders)
p

En las dos primeras se llama directamente a ggplot, bien definiendo los datos como argumento de la función, bien usando una tubería que «canalice» los datos desde el objeto murders hasta la función.

Las dos últimas líneas asocian el resultado de la función ggplot a una variable p. Para ver el gráfico (repito: aún no se ve nada) se muestra el valor de p11. Si no se asocia el gráfico a una variable, la función muestra el gráfico directamente, que es lo que sucede con cualquiera de las dos primeras líneas del ejemplo.

Una vez definido el objeto, los diferentes componentes se añaden usando capas. Las capas pueden «contener» la geometría, el aspecto, la escala de los ejes, calcular resúmenes estadísticos, etc. Para añadir capas se usa el símbolo «+», de manera que una línea de código de ggplot tiene un aspecto similar a data %>% ggplot() + capa1 + capa2 + capa3 + … + capan.

Normalmente la primera capa define la geometría, la que decide qué tipo de gráfico se va a representar. En el ejemplo necesitamos un gráfico de dispersión (scatterplot). Para conseguirlo, un vistazo a la chuleta de la función ggplot dice que la función necesaria es geom_point()12. Esta función necesita datos, que por defecto toma de los que hemos encaminado con la función «tubería» (pipe: %>%) a ggplot, y un mapeo o esquema. Este mapeo se puede consultar en la sección de estética (aesthetics) del archivo de ayuda: ?geom_point o help(geom_point). Los argumentos obligatorios son x e y, que son las coordenadas x e y del gráfico, ambos ejes.

Una forma de conectar los datos con el aspecto visual que se consigue en el gráfico es con la función aes (de aesthetics 'estética'). La salida de esta función se puede (y se suele) usar como elemento de entrada de la función geométrica (en este caso, geompoint).

El siguiente código proporciona un gráfico de dispersión en el que se ve el total de asesinatos del conjunto de datos del ejemplo comparado con el total de habitantes en millones:

murders %>% ggplot() +
  geom_point(aes(x = population/10^6, y = total))

El mismo código pero sin la función aes dará un error, en parte porque geompoint no entiende los argumentos population y total, pues son parte de un conjunto de datos, no variables definidas en el espacio de trabajo. Sustituyendo population por murders$population y total por murders$total desaparece el error, pero el gráfico sale en blanco. Por lo que se ve, geompoint() tiene como argumento la salida de lo que realmente dibuja el gráfico, que es aes. geompoint lo que hace es darle el aspecto deseado (gráfico de dispersión), que podría ser otro (geombars genera un gráfico de barras para la misma definición de datos, por ejemplo).

Otra cosa: en el código anterior se podrían obviar las etiquetas x = e y =, pues son los dos primeros argumentos, son obligatorios y están en el orden requerido.

Por último, las capas se pueden añadir a objetos que tengan datos asociados a ggplot. Por ejemplo, anteriormente se asignó la función ggplot a p y, por lo tanto, se puede añadir la capa a p así:

p + geom_point(aes(x = population/10^6, y = total))

En este caso, como se mencionó antes, no hay asignación a un objeto sino que se añade una capa a un objeto y el gráfico es mostrado inmediatamente. Para asociarlo a un objeto hay que asignarlo de nuevo: p1 <- p + geom_point(aes(x = ...)).

De momento, cosas como las etiquetas de los ejes, la escala, etc., son asignadas a sus valores por defecto.

Otra capa añade las etiquetas (texto) a cada uno de los puntos del gráfico, porque de momento no se sabe a qué estado corresponde cada punto del gráfico. Son las funciones geomlabel y geomtext, que también actúan sobre el mapeado estético de aes con sus correspondientes x e y obligatorias. El texto se añade con el argumento label en la función aes pasada a geomtext. En este caso se usa geomtext porque geomlabel hace lo mismo pero poniendo el texto en una cajita. Así que la cosa va quedando tal que así:

p + geom_point(aes(x = population/10^6, y = total)) +
geom_text(aes(x = population/10^6, y = total, label = abb))

Otra vez, la etiqueta label = abb no puede sacarse «fuera» de aes porque abb no es una variable definida en el espacio de trabajo, aunque podría hacerse escribiendo el nombre completo de la variable del conjunto de datos: murders$abb.

Aquí veo claramente la necesidad de no escribir todo el tiempo los valores de x e y, por más que sea posible obviar lo de x=, y escribir una variable o un comodín. Otro problema que aparece: el texto está escrito exactamente encima del punto correspondiente, con lo que el gráfico es un caos indescifrable.

El segundo problema se corrige con facilidad usando uno de los argumentos opcionales de geom_text(): la opción nudge_x, de nudge 'codazo, empujón', «aparta» un poco las etiquetas de los puntos en el eje x, así que el código es, ahora:

p + geom_point(aes(x = population/10^6, y = total), size = 3) +
geom_text(aes(x = population/10^6, y = total, label = abb), nudge_x = 1)

Nótense dos cosas: una, que en la función geom_point() está incluido el parámetro size = 3 para hacer los puntos algo mayores y más visibles; y dos, que tanto este parámetro como nudge están «fuera» de aes, no hay un mapeo entre los datos y el gráfico, sino simplemente se define un valor u opción.

El primer problema se soluciona del modo siguiente: definiendo un mapeado estético global cuando se asocia el conjunto de datos a ggplot13, dentro del argumento (los paréntesis) que se habían quedado en blanco. El nuevo código es, por lo tanto:

murders %>% ggplot(aes(population/10^6, total, label = abb))
p + geom_point(size = 3) +
geom_text(nudge = 1)

Los argumentos que controlan el tamaño de los puntos y la separación del texto continúan estando en sus respectivos lugares, ya que son específicos de esas funciones y no tienen sentido dentro de ggplot. Pero hay algo muy importante: los mapeados estéticos que se incluyan en las capas que se añaden sobrescriben los mapeados definidos globalmente, cosa ciertamente útil si se quiere modificar algo en una capa determinada. Por ejemplo, añadir un texto en la capa de etiquetas: + geom_text(aes(x = 10, y = 800, label = "Hola, mundo")). Ahora aparece un «Hola, mundo» centrado en el punto (10, 800) del gráfico en lugar de las etiquetas abb definidas en ggplot.

  1. Ajustar escalas, colores y textos

    Ya está listo el gráfico completo, pero estéticamente necesita muchos ajustes. El apelotonamiento de puntos cerca del centro de coordenadas hace que sea muy difícil de interpretar, por ejemplo.

    Una de las primeras cosas que hay que hacer es ajustar los ejes a una escala logarítmica, es decir, que en lugar de que los intervalos regulares sean los puntos 0, 10, 20, 30, etc., que sean los puntos 0, 10, 100, 1000… Para eso, existen las funciones scale_x_continuous() y scale_y_continuous(), una para cada eje porque ambos han de ser modificados por separado (se puede modificar uno y el otro no, si hace falta). Se añaden, por lo tanto, dos nuevas capas:

    scale_x_continuous(trans = "log10") + 
    scale_y_continuous(trans = "log10")
    

    Sin embargo, la transformación logarítmica es tan común que hay funciones específicas para ello (scale_x_log10()). El nuevo código será:

    murders %>% ggplot(aes(population/10^6, total, label = abb))
    p + geom_point(size = 3) +
    geom_text(nudge = 0.075) +
    scale_x_log10()
    scale_y_log10()
    

    Al hacer la escala logarítmica, el «empujoncito» que le dimos a las etiquetas abb es excesivo. Será mejor ajustarlo a 0.075, por ejemplo.

    Los ejes aún no tienen ningún texto que haga referencia a los datos que se están representando. Las funciones que se usan para los ejes y el título del gráfico son (incluyendo los valores válidos para el ejemplo):

    xlab("Población en millones (escala logarítmica)") + 
    ylab("Número total de asesinatos (escala logarítmica)")
    ggtitle("Asesinatos con arma de fuego en EEUU (2010)")
    

    En realidad, el gráfico ya se ve bastante bien y es comprensible e interpretable. Pero aún hay cosas que se pueden mejorar. Una es añadir colores a los puntos para codificar información como la región a la que pertenece el estado. El parámetro col dentro de geom_point() permite hacer esto, y se hará redefiniendo el objeto p sin los puntos, que se añadirán luego.

    murders %>% ggplot(aes(population/10^6, total, label = abb)) +
    geom_text(nudge = 0.075) +
    scale_x_log10() +
    scale_y_log10() +
    xlab("Población en millones (escala logarítmica)") + 
    ylab("Número total de asesinatos (escala logarítmica)") +
    ggtitle("Asesinatos con arma de fuego en EEUU (2010)")
    

    Con este p ya definido, se puede añadir color (azul) a los puntos con p + geom_point(size = 3, col = "blue"). Pero esto no aporta información, solo cambia el color. Para aportar información se mapean los datos de la variable region en la función geom_point():

    p + geom_point(aes(col = region), size = 3)
    

    El asociar colores a los valores de una variable categórica y mostrar una leyenda con esa asociación es el comportamiento por defecto de ggplot.

    Falta un tema importante. Fijándose un poco, resulta que la leyenda que muestra los colores de los puntos y su significado proviene de la variable regions y, como tal, tiene ese nombre. Pero no es un nombre apropiado, debería, al menos, comenzar con mayúscula. Una nueva línea de código permite personalizar el nombre de la leyenda:

    + scale_color_discrete(name = "Region")
    
  2. Línea media

    Una línea discontinua gris marcará dónde está la media de asesinatos en todo el país. Para calcular cuál es esa media se multiplica el ratio de asesinatos por cada millón de habitantes r por la población en millones. La fórmula será y = r*x.

    El código usado para este cálculo necesita la librería dplyr, que ya debería estar cargada para usar el operador %>%.

    r <- murders %>%
      summarize(rate = sum(total) / sum(population) * 10^6) %>%
      .$rate
    

    Explicación del código: se dirige el conjunto de datos a la función summarize, que suma el total de asesinatos y la población, divide ambos y multiplica el resultado por 106. Esto se deriva a una nueva variable rate en el conjunto de datos. Todo esto se asocia al nuevo objeto r.

    La línea se dibuja usando la función geomabline. Esta función «intercepta» el valor a con pendiente b. Los valores predeterminados son 0 y 1, respectivamente. En este caso solo es necesario explicitar el valor de intercept, que será el logaritmo decimal de r. Esto necesita más elaboración y una mejor explicación. En cualquier caso, el código para la línea será (es otra capa que se añade):

    + geom_abline(intercept = log10(r), lty = 2, color = "darkgrey")
    

    Directamente se especifica, también, que la línea sea discontinua con lty = 2 (de line type) y el color. Hay que tener en cuenta que es mejor añadir la capa de la línea antes de los puntos y el texto. De esta manera, la línea no se superpone y tanto los puntos como las etiquetas abb son mucho más legibles.

3.4 Resumir con dplyr

Dentro del paquete tidyverse está la función summarize, que ya apareció un poco más arriba (pero sin explicar nada). Esta función ofrece resúmenes de un conjunto de datos en una tabla resumen. Los nombres de los datos resumidos son definidos en la propia función.

Como la mayoría de las funciones del paquete dplyr/tidyverse, la función summarize es capaz de entender los nombres de las variables dentro del conjunto de datos. Es decir, que no es necesario escribir el conjunto de datos, el accesor $ y el nombre de la variable.

El resultado de summarize es en sí un pequeño conjunto de datos a cuyas variables se puede acceder igual que a las de cualquier otro conjunto de datos. El siguiente código hace lo siguiente: pasa el conjunto de datos heights como parámetro a la función filter para extraer solo los valores correspondientes a la variable sex cuando es igual a Male (es decir, solo los varones), y, a su vez, pasa este resultado como parámetro a la función summarize para elaborar una tabla resumen, la cual asigna al objeto s. A continuación muestra las variables dentro de s que han sido definidas en la misma función summarize.

s <- heights %>%
    filter(sex == "Male") %>%
    summarize(media = mean(height), sigma = sd(height))

# acceder a la media (media) y desviación típica (sigma) desde la
# tabla resumen
s$media
s$sigma

En un principio, summarize no puede usarse con funciones que devuelven varios valores, como quantile, sino solamente con funciones que devuelven un valor único, como median, min, sd, etc. Sin embargo, a partir de la versión 1.0 de dplyr esto ya no es así. Por ejemplo, summarize(rango = quantile(height, c(0, 0.5, 1))) no devolverá un error si dplyr es 1.0 o superior.

3.4.1 El marcador punto

En ocasiones, puede ser interesante que las funciones de dplyr devuelvan vectores en lugar de data frames. Un ejemplo sería usar la función summarize para calcular la media nacional de asesinatos por cada 100,000 habitantes. Recapitulando: hemos importado los datos y añadido una columna con el ratio por 100,000 habitantes. Este es el conjunto de datos más completo que tenemos hasta ahora:

data(murders)
murders <- murders %>% mutate(murder_rate = total/population*100000)

A continuación, calculamos la media de los 50 valores de la variable murderrate:

summarize(murders, mean(murder_rate))

Lo cual devuelve el valor 2.779125, que es incorrecto. ¿Por qué? Porque no tiene en cuenta el valor de la población de cada estado, tan solo es una media de los ratios. Los estados más pequeños tienen el mismo peso que los más grandes y poblados. La media nacional debe calcularse sumando las poblaciones y sumando los totales de asesinatos, no operando directamente con los ratios, así:

us_murder_rate <- murders %>%
  summarize(rate = sum(total) / sum(population) * 100000)

Ahora sí, el resultado es 3.03455.

Esto, en realidad, es solo una digresión para mostrar qué es exactamente la variable murderrate y cómo es manejada por R. El problema viene ahora: el resultado de summarize es un número (3.03455), pero class(us_murder_rate) devolverá data frame. Entonces, si necesitamos operar con ese 3.03455, ¿qué podemos hacer?

Pues acceder al valor en sí dentro de la tabla, y para eso se puede escribir tanto el nombre completo us_murder_rate$rate como usar el marcador u operador punto. usmurderrate %>% .$rate` es equivalente al nombre completo. El punto funciona exactamente igual que el operador del directorio actual en la consola de Bash14. Ejemplo:

r <- us_murder_rate %>% .$rate
r
[1] 3.03455

class(r)
[1] "numeric"

Este código devuelve el valor de la variable rate dentro del conjunto de datos usmurdersrate, que fue definido más arriba pasando murders a la función summarize, en la que se define una variable rate como el resultado de la operación de total entre population multiplicada por 10,000. Por lo tanto, se puede hacer lo mismo en una sola línea de código:

us_murder_rate <- murders %>%
    summarize(rate = sum(total) / sum(population) * 100000) %>%
    .$rate

Ahora, class(us_murder_rate) ya es numeric, pues el valor de la variable se ha pasado al objeto. Es un solo valor, pero podían ser varios, lo importante es que ahora es un vector con el que se puede operar numéricamente y no un conjunto de datos.

Ojo con que la variable detrás de $ es la misma que se define en summarize. Es una nueva variable en el conjunto de datos que acaba de ser definida y es la que se está convirtiendo de tabla en vector.

3.4.2 Agrupar y resumir

Una operación sumamente común es la de separar los datos en grupos y luego elaborar una visión resumida para cada grupo. Por ejemplo, calcular la media de altura para hombres y mujeres por separado.

Esto se puede lograr con la función groupby.

Si se canaliza el conjunto de datos heights en la función groupby y se especifica la variable sex como parámetro de agrupación, se obtiene una tabla. Esta tabla es un tipo especial denominado group data frame o conjunto de datos agrupados. Las funciones de dplyr, como summarize, se comportarán diferente en este tipo de tablas.

Pero cuando se ejecuta summarize, ofrece los datos de ambos grupos separadamente (ambos grupos porque en este conjunto de datos son dos grupos: male y female; en otros pueden ser más - por ejemplo, en data(murders), la orden group_by(region) ofrecerá 4 grupos).

heights %>% group_by(sex) %>%
  summarize(media = mean(height), sigma = sd(height))

Este código mostrará un resumen (como los anteriores) pero con dos filas en lugar de una, una para cada grupo.

Ejemplo con el conjunto de datos murders:

murders <- mutate(murders, ratio = total/population*100000)
murders %>% group_by(region) %>%
+ summarize(mediana = median(ratio))

Ahora la tabla muestra cuatro filas, una por región, con la mediana del ratio de asesinatos por cada 100,000 habitantes para cada región.

Nótese que groupby permite agrupar más de una variable, de manera que group_by(variable_1, variable_2) funcionará perfectamente, agrupando los datos en esas dos categorías.

3.4.3 Ordenar tablas (arrange)

Ordenar el contenido de una tabla para extraer algún conocimiento acerca de los datos contenidos es muy útil. Aunque ya se vieron las funciones sort y order para ordenar vectores, la función para ordenar tablas es arrange. Está contenida en dplyr.

Por ejemplo, el siguiente código ordena el conjunto de datos murders por población del estado: murders %>% arrange(population) %>% head(). Se toma el conjunto de datos y se canaliza hacia arrange, donde se especifica la columna que va a ser ordenada. El conjunto de datos completo se ordena en función de esa columna. La salida de esta función se canaliza a su vez hacia head, para mostrar las primeras 6 filas y sus encabezados nada más en lugar de las 51 entradas del conjunto de datos. Para ver más datos (pero no todos) se puede usar la función topn. Por ejemplo, top_n(10) muestra las primeras 10 entradas del conjunto de datos: murders %>% arrange(population) %>% top_n(10). La función topn15 puede seleccionar los valores máximos por sí sola, pero no ordenarlos: murders %>% top_n(10, ratio) mostrará los 10 primeros estados cuyo ratio es mayor (pero sin ordenar).

El comportamiento predeterminado es ordenar la tabla por la variable escogida siempre de menor valor a mayor (orden ascendente), pero se puede modificar usando la función desc (orden descendente), como en murders %>% arrange(desc(ratio)) %>% head().

Otra manera de ordenar una muestra del conjunto de datos es usar slicemax o slicemin, que escogen una muestra configurable de valores ordenados de mayor a menor o al revés. Así que slice_max(murders, n = 10, order_by = ratio) ofrece las primeras 10 entradas del conjunto de datos que tengan mayor ratio de asesinatos (la variable ratio usada para el ejemplo ya se definió en un punto anterior, creo que cuando se explicaba la función mutate), ordenadas de mayor a menor. De igual manera, slice_min(murders, n = 5, order_by = population) muestra una tabla con todos los datos de los estados con menor población, ordenados de menos a más.

  1. Ordenamiento anidado

    También se pueden ordenar tablas anidando criterios en lugar de usar un criterio simple. Por ejemplo, en el conjunto de datos murders se puede ordenar por región y luego, dentro de la región, por ratio de asesinatos por cada 100,000 habitantes. Simplemente, se especifican los criterios en orden (por supuesto, esto es compatible con la orden desc) igual que antes:

    murders %>% arrange(region, ratio) %>% head()
    

    En estos criterios de ordenación anidados se puede, como se dijo antes, usar desc. Pero desc atañe solo al primer parámetro, es decir, que arrange(desc(region, population)) ordena las regiones pero no ordena los estados por población dentro de las regiones. Para eso es necesario añadir otro desc: arrange(desc(region), desc(population)). Así, el orden de los estados dentro de cada region estará ordenado de mayor a menor.

    Un ejemplo algo diferente para probar un conjunto de datos distinto en el que se pueden usar varios criterios. El objetivo es conocer las presiones arteriales medias y las desviaciones típicas, por raza, ordenadas de menor a mayor, de solamente los varones entre 40 y 49 años. El conjunto de datos hay que instalarlo e importarlo, y se refiere a estudios realizados en los EEUU.

    install.packages("NHANES")
    library(NHANES)
    library(dplyr)
    data(NHANES)
    
    NHANES %>%
      filter(Gender == "male" & AgeDecade == " 40-49") %>%
      group_by(Race1) %>%
      summarize(media = mean(BPSysAve, na.rm = TRUE),
        sigma = sd(BPSysAve, na.rm = TRUE)) %>%
      arrange(media)
    

3.5 Gapminder

3.5.1 El conjunto de datos Gapminder

Esta sección es un caso práctico. «Gapminder» es un conjunto de datos incluido en el paquete dslabs, y contiene datos de una hoja de cálculo disponible en el sitio gapminder.org. Es un sitio de internet que utiliza datos reales y buenas técnicas de visualización para desmontar mitos y falacias a propósito del estado actual de las cosas en cuanto a la pobreza en el mundo, el cambio climático, etc.

Como se dijo, está incluido en dslabs, de manera que para poder usar este conjunto de datos:

library(dslabs)
data(gapminder)

Usando estos datos se pueden extraer conclusiones realistas de un simple vistazo, siempre y cuando los gráficos estén bien definidos.

Un ejemplo es el uso de la tasa de fertilidad (el número de hijos por mujer) y la esperanza de vida. La idea preconcebida es que el mundo se puede dividir, grosso modo, en dos: el mundo desarrollado, con baja natalidad y alta esperanza de vida (Europa Occidental, Norteamérica y Japón) y un mundo en desarrollo (mayoritariamente los países de África, Asia y Sudamérica). Para comprobar si esta visión dicotómica del mundo se corresponde o no con la realidad se hará un gráfico con la situación de hace 50 años y luego uno con la situación actual:

library(dslabs)
data(gapminder)
library(dplyr)
library(ggplot2)
ds_theme_set()
filter(gapminder, year == 1962) %>% 
  ggplot(aes(fertility, life_expectancy)) +
  geom_point()

Los puntos se distribuyen en dos áreas más o menos diferenciadas, pero si se añade un código de color para el continente en el que se encuentra cada uno de los puntos, se ve claramente que, efectivamente, Europa está separada del resto. Para añadir el código de colores: ggplot(aes(fertility, life_expectancy, color = continent)) + geom_point(). El gráfico resultante es:

gapminder_1.png

Figure 5: Distribución de países por continente en 1962.

Para saber si esta distribución polarizada sigue siendo real en la actualidad se puede repetir el gráfico cambiando el año en el código, de 1962 a 2012, último año registrado en el conjunto de datos. Un vistazo al nuevo gráfico pone de manifiesto alguna dificultad:

gapminder_2.png

Figure 6: Distribución de países por continente en 2015

La más importante (de las dificultades) es que los ejes no tienen un rango coherente, por lo que hay que fijarse bien en las cifras para no equivocarse y poder interpretar algo. Otra dificultad es que hay que andar cambiando de uno a otro constantemente para comparar.

Para poder comparar gráficos como está mandado y apreciar diferencias de un vistazo hay que aprender acerca del «tallado» o faceting, es decir, elaborar varios gráficos (2 o más) que estén juntos, alineados y en una estructura coherente, como las caras de un diamante u otra gema tallada (de ahí el nombre, supongo16).

3.5.2 Faceting o tallado

Hay varias funciones que permiten organizar varios gráficos de manera que sea fácil compararlos. Dependiendo de las necesidades del momento se puede usar una u otra.

El faceting facilita la comparación de datos segmentando los datos por una variable determinada y haciendo el mismo gráfico varias veces para cada valor de esa variable. En el caso del ejemplo, perteneciente al conjunto de datos Gapminder, se puede hacer esto año a año.

La función que permite el tallado o faceting es facet_grid(). Esta función permite comparar hasta dos variables a la vez usando las filas para una variable y las columnas para la otra, y se añade como una capa más al código del gráfico. Las filas y columnas se separan con una tilde (~): facet_grid(continent~year):

filter(gapminder, year%in%c(1962, 2012)) %>%
  ggplot(aes(fertility, life_expectancy, col = continent)) +
  geom_point() +
  facet_grid(continent~year)

En este ejemplo, esta información puede ser excesiva. Podemos querer comparar tan solo dos años, en dos columnas, para poder ver los gráficos uno al lado del otro. En tal caso, solo se especifica una variable en lugar de dos, y el punto se usa para sustituir la variable que no se necesita, de manera que todavía es posible usar la tilde para definir filas y columnas.

filter(gapminder, year%in%c(1962, 2012)) %>%
  ggplot(aes(fertility, life_expectancy, col = continent)) +
  geom_point() +
  facet_grid(.~year)

Lo cual da como resultado:

facet-0.png

Figure 7: Faceting de una variable

Sin embargo, otro problema aparece en el horizonte. Dado lo interesante del gráfico, sería una buena idea añadir más años a la lista para ver la evolución. Pero en tal caso, colocar todos los gráficos en una sola fila los hará demasiado estrechos para ser útiles. Para eso se usa la función facetwrap.

años <- c(1962, 1980, 1990, 2000, 2012)
filter(gapminder, year%in% years) %>%
  ggphot(aes(fertility, life_expectancy, col = continent)) +
  geom_point() +
  facet_wrap(~year)

Ahora es más fácil ver la evolución temporal de la tasa de fertilidad y la esperanza de vida.

Es muy importante recordar que los ejes de los gráficos no están definidos por los datos en el gráfico, sino por los datos en todos los gráficos, de manera que es sencillo comparar al primer golpe de vista. Si los gráficos se dibujasen uno a uno, los ejes no mostrarían el mismo rango y sería más complicado comparar.

3.5.3 Series temporales

Los gráficos de puntos pueden mostrar la evolución temporal de una variable (es decir, los valores de la variable en función de otra variable de tipo tiempo, como years en este conjunto de datos del ejemplo). Sin embargo, es mucho más conveniente mostrar esta evolución como una línea, que no es más que la unión de los puntos.

Para dibujar líneas solo hay que cambiar la capa geompoint por la capa geomline.

Esto es particularmente útil cuando se comparan dos o más variables (por ejemplo, dos o más países o regiones en Gapminder). Pero no funciona así, tal cual. Si se pasan a ggplot dos países (o más), el gráfico intenta unir todos los puntos, independientemente de cuál sea el valor o país que representa y el resultado es un desastre. Hay que mapear la variable que tiene los valores que se desean comparar (por ejemplo, si son dos o más países -Alemania contra Corea-, la variable country). Un código de ejemplo con Alemania y Corea, ya que estamos:

paises <- c("Germany", "South Korea")
gapminder %>% filter(country %in% paises) %>%
  ggplot(aes(year, fertility, group = country)) +
  geom_line()

Una posible mejora es sustituir el parámetro group dentro del mapeo estético por color o col. Group efectivamente agrupa los datos en dos líneas separadas, pero no informa nada acerca de qué línea corresponde a cada país, pero si se incluye el parámetro color, así: ggplot(aes(year, fertility, color = country)) se matan dos pájaros de un tiro (en realidad, tres: datos agrupados, líneas coloreadas y una leyenda explicativa).

Otra mejora es sustituir la leyenda por defecto con etiquetas dentro del mismo gráfico. Este es el modo recomendado de representar los datos. Para hacer eso hay que definir un conjunto de datos con las etiquetas y usar un mapeo de texto para colocar las etiquetas en su sitio. El ejemplo a continuación usa la variable lifeexpectancy 'esperanza de vida':

# El objeto 'paises' ya ha sido definido antes
etiquetas <- data.frame(country = paises, x = c(1975, 1965), y =
c(72, 60))
gapminder %>% filter(country == paises) %>%
  ggplot(aes(year, life_expectancy, col = country)) +
  geom_line() + 
  geom_text(data = etiquetas, aes(x, y, label = country), size = 5) +
  theme(legend.position = "none")

Una cosa importante a tener en cuenta aquí es que los valores de x e y que se definen en el objeto etiquetas están puestos a ojo. Posiblemente haya una manera de hacerlo automáticamente, pero en principio hay que ir tanteando.

3.5.4 Transformaciones: entender mejor las distribuciones

Las transformaciones pueden ser muy útiles para entender mejor las distribuciones. En este apartado se usarán los datos de ingresos económicos disponibles dentro del conjunto de datos Gapminder en la variable gdp, Gross Domestic Product por sus siglas en inglés o PIB en román paladino. El PIB per cápita es una estimación usual aproximada de la riqueza de un país. Si se divide por 365 se obtiene el ingreso diario per cápita. Este dato se puede interpretar de la siguiente manera:

  • Menos de $2: pobreza extrema
  • 2$: muy pobre
  • 4$: pobre
  • 8$: riqueza media
  • 16$: país con recursos
  • 32$: país rico
  • 64$: muy rico

Por lo tanto, se incluirá una nueva variable en el conjunto de datos con este valor:

gapminder <- gapminder %>%
  mutate(dolares_diarios = gdp/population/365)

Ahora hay una nueva columna llamada dolaresdiarios que indica la renta per cápita diaria en dólares americanos actuales. Un sencillo histograma muestra los valores de la variable dólaresdiarios:

año_pasado <- 1970
gapminder %>% filter(year == año_pasado & !is.na(gdp)) %>%
  ggplot(aes(dolares_diarios)) +
  geom_histogram(binwidth = 1, color = "black")

La mayoría de países está por debajo de 10 dólares diarios, aunque la mayor parte del eje X está dedicado a los países por encima de 10 $, que son minoría. Dado que los valores que consideramos útiles son potencias de 2 (1$ diario, 2$ diarios, 4$ diarios, 8$ diarios, etc.), una transformación de los datos para mostrar el logaritmo en base 2 en lugar del valor real ofrecerá un histograma completamente diferente: ggplot(aes(log2(dolares_diarios))) (el resto del código es igual). El histograma luce tal que así:

log2-histo.png

Figure 8: Histograma del logaritmo en base 2

Se pueden observar dos picos en lugar de uno. Cada uno de los picos se denomina moda. La moda es el valor que se repite con más frecuencia, y en una distribución normal coincide con la media. Cuando los valores no decrecen suavemente se habla de «modas locales», y aquí hay dos: una alrededor de 2$ diarios y otra alrededor de 32$ diarios (coincide con el valor 5 en el eje X, es decir, 25 o 32). Esta distribución bimodal en 1970 centrada en 2$ y 32$ es consistente con la separación del mundo en dos ejes: países pobres y países ricos.

La elección de los logaritmos en base 2 o 10 son usuales y muy útiles, pero no se recomienda en absoluto el uso de logaritmos naturales (potencias de e) o neperianos, pues es muy difícil interpretar con un cálculo mental los valores reales representados en el eje X. A la hora de elegir la escala, es mejor y más fácil cuantos más números enteros haya y menos números decimales. En el ejemplo la escala está entre -2.5 y 6.5, pero si se hubiera escogido el logaritmo en base 10, los valores estarían entre 0 y 1 y es más complicado interpretar potencias expresadas de esta manera. La elección del ancho de las barras también es más sencilla. Un ejemplo en el que el logaritmo en base 10 es más adecuado sería para representar la población, pues el rango está entre 45 000 y 800 millones (aprox.).

La transformación logarítmica se puede hacer de dos maneras. Ambos enfoques son útiles y tienen sus propias ventajas:

  1. Transformar los datos y representarlos: es más fácil interpretar los ejes, pues la escala es natural y continua, de manera que los valores intermedios se pueden deducir con suma facilidad (el punto medio entre 0 y 5 es 2.5 y se ve a simple vista).
  2. Representar los datos en un eje con una escala logarítmica: los valores en el eje X son más difíciles de interpretar por lo contrario en el punto anterior. En una escala 0 - 10 - 100 - 1000, el punto medio entre 10 y 100 no es 50. Como ventaja, los valores en el eje X son valores reales, no potencias de 2 o de 10 que haya que calcular para leer el gráfico.

Para implementar una escala logarítmica se usan los ajustes de los ejes ya vistos anteriormente: scale_x_continuous y se dejan los datos tal cual, no como en el ejemplo anterior. El mismo histograma con un eje logarítmico sería:

año_pasado <- 1970
gapminder %>% filter(year == año_pasado & !is.na(gdp)) %>%
  ggplot(aes(dolares_diarios)) +
  geom_histogram(binwidth = 1, color = "black") +
  scale_x_continuous(trans = "log2")

Donde antes había un 5 en el eje X ahora hay un 32, porque 25=32, es decir, 32$ dólares diarios. Ahora ya se sabe, pero a primera vista no es tan evidente que cada barra representa los valores 1, 2, 4, 8, 16, 32 y 64, por lo que leer los valores de barras sin un número debajo se complica un poco.

3.5.5 Estratificación y diagramas de cajas y bigotes

Los ejemplos anteriores muestran una división del mundo en 1970 en dos polos. Pero no dicen nada acerca de si es una división Este-Oeste, países ricos-países en desarrollo, etc. Habría que hacer un histograma por cada región en que se divide el conjunto de datos y comparar, pero dado que hay 22 regiones diferentes, se hace un poco difícil y poco práctico.

En su lugar, se usarán diagramas de cajas y bigotes unos junto a otros:

p <- gapminder %>% filter(year == año_pasado & !is.na(gdp)) %>%
  ggplot(aes(region, dolares_diarios))
p + geom_boxplot()

La información está toda ahí, pero hay varios problemas con el gráfico. El más evidente es que los valores del eje X (region) se superponen unos a otros. Para rotar las etiquetas y que sean legibles se modifica la última línea del código añadiendo otra capa: theme(axis.text.x = element_text(angle = 90, hjust = 1))17.

Hay otros cambios necesarios: ordenar los datos, agregar colores para codificar los continentes a los que pertenecen los países y ajustar la escala del eje Y a una escala logarítmica. Para reordenar se usa la función reorder asociada a la variable region dentro de la función mutate: se redefine la función region reordenándola basándose en la variable dolaresdiarios y usando la función median.

p <- gapminder %>$
  filter(year == año_pasado & !is.na(gdp)) %>%
  mutate(region = reorder(region, dolares_diarios, FUN = median)) %>%
  ggplot(aes(region, dolares_diarios, fill = continent)) +
  geom_boxplot() +
  scale_y_continuous(trans = "log2") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

En casos como este, en los que el número de datos concretos no es muy alto, se pueden incluir los puntos superpuestos al gráfico añadiendo una capa geompoint eliminando las etiquetas: geom_point(show.legend = FALSE). En cualquier caso, hay que ser muy cuidadoso y no sobrecargar la imagen. Demasiada información puede hacer que todo sea más difícil de leer y se pierda capacidad de síntesis.

3.5.6 Comparar distribuciones

Hasta ahora se ha visto que la distribución de la riqueza es una distribución bimodal (en el histograma) y, mirando los diagramas de cajas y bigotes anteriores, que los países ricos están en Occidente en oposición al resto del mundo. Para corroborar esta hipótesis, hay que hacer histogramas en los que se separen los países occidentales del resto.

Para ello se crea un vector que incluya solo los países del mundo occidental: occidente <- c("Western Europe", "Northern America", "Southern Europe", "Northern Europe", "Australia and New Zealand"). Después se filtran los países siguiendo este vector y se generan los gráficos:

gapminder %>%
    filter(year == año_pasado & !is.na(gdp)) %>%
    mutate(grupo = ifelse(region %in% occidente, "Occidente", "En desarrollo")) %>%
    ggplot(aes(dolares_diarios)) +
    geom_histogram(binwidth = 1, color = "black") +
    scale_x_continuous(trans = "log2") +
    facet_grid(. ~ grupo)

Lo interesante de este código está en el uso de una condición ifelse anidada en una función mutate. Lo que hace es decidir si el valor de region está contenida en el vector occidente; si es así, el valor de la variable grupo es Occidente; si no, es En desarrollo; y así completa todas las entradas de la tabla. Después se genera un histograma de la renta per cápita y se establece una variable (grupo) para crear una «malla» de histogramas (la variable grupo definida en la segunda línea se muestra en columnas, es decir, los histogramas se colocan uno al lado del otro) usando facetgrid.

histo1.png

Figure 9: Histogramas para comparar Occidente al resto del mundo

Está claro que los países occidentales son más ricos que el resto del mundo, pero ¿será cierto lo mismo en 2010? Se hace lo mismo pero incluyendo los años a estudiar en un vector que incluya ambos objetos previamente definidos (dentro de la función filter) y se incluye la variable year dentro de la parrilla de histogramas a generar:

año_actual <- 2010
gapminder %>%
    filter(year %in% c(año_pasado, año_actual) & !is.na(gdp)) %>%
    mutate(grupo = ifelse(region %in% occidente, "Occidente", "En desarrollo")) %>%
    ggplot(aes(dolares_diarios)) +
    geom_histogram(binwidth = 1, color = "black") +
    scale_x_continuous(trans = "log2") +
    facet_grid(year ~ grupo)

dos puntos del tiempo histo2.png

Nuevo problema: prestando un poco de atención se ve que las barras muestran un conteo de países superior en 2010 que en 1970. Lógico, hay países nuevos, como los creados tras la caída de la Unión Soviética en los 90 y hay más datos ahora que en 1970. Lo debido sería comparar solo los países de los cuales se tienen datos en ambos puntos del tiempo:

lista_1970 <- gapminder %>%
  filter(year == año_pasado, !is.na(gdp)) %>%
  .$country

lista_2010 <- gapminder %>%
  filter(year == año_actual, !is.na(gdp)) %>%
  .$country

lista_paises <- intersect(lista_1970, lista_2010)

Este código crea una lista en la intersección entre dos listas previas. Estas listas se crean filtrando los datos por año (como hasta ahora), y pasando el vector lógico que se crea a la variable country, de manera que se tiene una lista de países para los que sí hay datos en 1970. Lo mismo para 2010. Y finalmente se usa la función intersect para unir ambas listas. Solo los países que pertenezcan a ambas listas formarán parte de la lista final.

Se rehacen los histogramas.

gapminder %>%
    filter(country %in% lista_paises & year %in% c(año_pasado, año_actual) &
	   !is.na(gdp)) %>%
    mutate(grupo = ifelse(region %in% occidente, "Occidente", "En desarrollo")) %>%
    ggplot(aes(dolares_diarios)) +
    geom_histogram(binwidth = 1, color = "black") +
    scale_x_continuous(trans = "log2") +
    facet_grid(year ~ grupo)

Está claro que los países occidentales han mejorado, pero los países en desarrollo han mejorado también (probablemente más, pero no está claro18).

Una manera estupenda de apreciar las diferencias entre ambos años sería poner las dos series de datos juntas en un diagrama de cajas y bigotes. Dado que ggplot colorea automáticamente las cajas dependiendo del valor en variables de tipo factor, lo que hay que hacer es indicarle a ggplot que year no es un número sino una variable de tipo factor.

gapminder %>%
    filter(year %in% c(año_pasado, año_actual) & !is.na(gdp)) %>%
    mutate(region = reorder(region, dolares_diarios, FUN = median)) %>%
    ggplot(aes(region, dolares_diarios, fill = factor(year))) +
    geom_boxplot() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    scale_y_continuous(trans = "log2")

boxpl0.png

Figure 10: Diagrama de cajas y bigotes por años

En lugar de hacer una tabla con varios gráficos, se incluyen todos los años pero se separan por colores. Así es más fácil comparar ambas épocas y el objetivo es hacer fácil la comparación.

3.5.7 Diagramas de densidad

Hasta aquí todo bien. Pero hay una manera de presentar toda esta información en un solo gráfico usando diagramas de densidad (algo de esto so vio antes durante la explicación de las distribuciones y los histogramas).

Viendo los histogramas está claro que la brecha entre países se cierra: pasamos de una distribución bimodal a una en la que esa diferencia se difumina. Para apoyar la tesis de que los países pobres se eriquecen en lugar de que haya países ricos que se empobrecen, hay que dar color a las distintas áreas en un gráfico de densidad.

Lo primero que hay que saber acerca de las curvas de densidad (ya se habló de esto antes) es que el eje Y no muestra la cuenta de elementos sino su frecuencia de aparición, haciendo que el área bajo la curva sea igual a 1 (es decir, su probabilidad de aparición). En este caso, dado que ambos grupos (Occidente y En desarrollo) tienen tamaños dispares (21 y 87, respectivamente) y que las áreas son de valor 1, las curvas no son equivalentes ni comparables.

Para que las curvas sean comparables, hay que multiplicar los valores del eje Y por el tamaño del grupo. La función geomdensity debe ser capaz de acceder a las variables y operar con ellas. Un vistazo a help(geom_density) muestra que hay una variable calculada (computed variable) que hace esto mismo: multiplicar la densidad por el número de puntos. Es la variable count. Para acceder a estas variables calculadas hay que rodearlas de la secuencia .., así: ..count...

gapminder %>%
    filter(country %in% lista_paises & year %in% c(año_pasado, año_actual) &
	  !is.na(gdp)) %>%
    mutate(grupo = ifelse(region %in% occidente, "Occidente", "En desarrollo")) %>%
    ggplot(aes(dolares_diarios, y = ..count.., fill = grupo)) +
    geom_density(alpha = 0.2) +
    scale_x_continuous(trans = "log2") +
    facet_grid(year ~ .)

Recuérdese que este código es solo el código del gráfico. Las librerías que hay que cargar y los objetos definidos anteriormente como añopasado o occidente se obvian porque en el mismo espacio de trabajo no hay que cargar o definir las cosas más de una vez. Pero si no están definidas, no funcionará.

Nueva vuelta de tuerca: en la exploración inicial de los datos se vio que, entre los países en desarrollo, son los asiáticos los que más han mejorado. ¿Cómo podemos mostrar esta información? Segmentando por regiones clave y coloreando diferentes curvas, una por región. Una nueva función llamada casewhen será muy útil para definir grupos. casewhen no acepta un argumento de datos, así que hay que mostrarle las variables directamente, con su nombre completo.

Vamos pues a definir los grupos:

gapminder <- gapminder %>%
  mutate(grupo = case_when(
    .$region %in% occidente ~ "Occidente",
    .$region %in%
      c("Eastern Asia", "South-Eastern Asia") ~ "Asia Oriental",
    .$region %in% c("Caribbean", "South America", "Central America") ~
      "Latinoamérica"
    .$continent == "Africa" &
      .$region != "Northern Africa" ~ "África subsahariana",
    TRUE ~ "Otros"
    ))

El código es autoexplicativo: cuando el contenido de region esté incluido en el vector occidente, el valor de grupo será Occidente, etc. Es interesante fijarse en que la lógica puede ser compleja: si el valor de continent es Africa y el valor de region no es Northern Africa, entonces se trata del África subsahariana.

Luego se convierte la variable grupo en una variable de tipo factor:

gapminder <- gapminder %>%
  mutate(grupo = factor(grupo, levels = c("Otros", "Latinoamérica",
    "Asia Oriental", "África subsahariana", "Occidente")))

Al generar el gráfico de nuevo, este es un poco difícil de leer. Para simplificar un poco las cosas, se usará la opción stack, que apila las curvas para que se vean mejor. Es una nueva opción dentro de la función geomdensity: geom_density(alpha = 0.2, bw = 0.75, position = "stack").

La última cosa necesaria para mostrar los datos correctamente es ponderar el tamaño de los países. Un país es un país, y computa lo mismo China o India, con más de 1,000 millones de habitantes, que países pequeños. Por lo tanto, hay que añadir otra opción weight al mapeado de ggplot. Antes, habrá que definir la operación necesaria para esa ponderación: population/sum(population)*2. No tengo nada claro por qué la ponderación se hace de esta manera; en fin, es la relación entre la población del país y la población total, pero… ¿y ese multiplicar por 2? Luego se incluye la opción weight = weight en el mapeado de ggplot. Como esto aún no está nada claro porque hay más modificaciones que no sé a qué vienen, se queda así. Ya veremos, el libro no da más explicaciones (apartado 9.7.5).

Al final, el código necesario para generar un gráfico como este desde cero y guardarlo en disco es el siguiente:

library(dslabs)
library(ggplot2)
library(dplyr)
data(gapminder)

# modificar el conjunto de datos con una nueva variable
gapminder  <- gapminder %>% mutate(dolares_diarios = gdp/population/365)

# definir un año a mostrar
año_pasado <- 1970
año_actual <- 2010

# Separar Occidente del resto del mundo
occidente  <- c("Western Europe", "Northern America", "Southern Europe", "Northern Europe", "Australia and New Zealand")


# guardar el gráfico generado en un archivo
png(file = "/home/javier/documentos/estudios_y_cursos/informatica/gnu_r/harvardx/apuntes/incl/grafico_de_prueba.png")

# definir las listas de países
lista_1970 <- gapminder %>%
  filter(year == año_pasado, !is.na(gdp)) %>%
  .$country

lista_2010 <- gapminder %>%
  filter(year == año_actual, !is.na(gdp)) %>%
  .$country

lista_paises <- intersect(lista_1970, lista_2010)

# agrupar los países en una nueva variable *grupo* y hacerla de tipo factor
gapminder <- gapminder %>%
  mutate(grupo = case_when(
    .$region %in% occidente ~ "Occidente",
    .$region %in%
      c("Eastern Asia", "South-Eastern Asia") ~ "Asia Oriental",
    .$region %in% c("Caribbean", "South America", "Central America") ~
	"Latinoamérica",
    .$continent == "Africa" &
      .$region != "Northern Africa" ~ "África subsahariana",
    TRUE ~ "Otros"
    ))

gapminder <- gapminder %>%
  mutate(grupo = factor(grupo, levels = c("Otros", "Latinoamérica",
    "Asia Oriental", "África subsahariana", "Occidente")))

# generar el gráfico
grafico  <- gapminder %>%
    filter(country %in% lista_paises & year %in% c(año_pasado, año_actual) &
	  !is.na(gdp)) %>%
    ggplot(aes(dolares_diarios, y = ..count.., fill = grupo)) +
    geom_density(alpha = 0.3, bw = 0.75, position = "stack") +
    scale_x_continuous(trans = "log2") +
    facet_grid(year ~ .)

# mostrar el gráfico
grafico

# devolver el control a la consola
dev.off()

3.5.8 La falacia ecológica

En toda la sección anterior se han estudiado las distribuciones de riqueza y otros parámetros en distintas regiones del mundo y su evolución en el tiempo. Pero los datos han de ser tomados con precaución para no caer en la https://es.wikipedia.org/wiki/Falaciaecol%C3%B3gica. Básicamente se trata de tener en cuenta que las medias por grupo no son representativas de los elementos tomados por separado. Se proponen a continuación un par de ejemplos para ilustrarlo usando, cómo no, unas cuantas funciones de R. El objetivo es investigar la relación entre ingresos y mortalidad infantil.

Primero, se definen unas cuantas regiones más a estudiar usando la función casewhen:

gapminder <- gapminder %>%
  mutate(grupo = case_when(
    .$region %in% occidente ~ "Occidente",
    .$region %in% "Northern Africa" ~ "Norte de África",
    .$region %in% c("Eastern Asia", "South-Eastern Asia") ~ "Asia Oriental",
    .$region == "Southern Asia" ~ "Sur de Asia",
    .$region %in% c("Caribbean", "South America", "Central America") ~
	"Latinoamérica",
    .$continent == "Africa" &
    .$region != "Northern Africa" ~ "África subsahariana",
    .$region %in% c("Melanesia", "Micronesia", "Polynesia") ~ "Islas del Pacífico"
    ))

Y luego se calculan los datos necesarios:

superv_ingresos  <- gapminder %>%
    filter(year %in% año_actual & !is.na(gdp) & !is.na(infant_mortality) &
	   !is.na(grupo)) %>%
    group_by(grupo) %>%
    summarize(ingresos = sum(gdp)/sum(population)/365,
	      ratio_superv_infantil = 1-sum(infant_mortality/1000*population)/sum(population))

La tasa de supervivencia infantil (a la que le hemos dado el nombre de ratiosupervinfantil) es lo contrario de la tasa de mortalidad, que es el dato que hay en el conjunto de datos. Por eso se calcula la mortalidad por 1,000 habitantes en relación a la población total y se resta de 1.

Una muestra de la tabla conseguida ordenada por la variable ingresos: superv_ingresos %>% arrange(ingresos). Es muy fácil ver en esta tabla la correlación entre ingresos y tasa de supervivencia infantil. Aún más claro si se dibuja el gráfico:

superv_ingresos %>%
    ggplot(aes(ingresos, ratio_superv_infantil, color = grupo, label = grupo)) +
    scale_x_continuous(trans = "log2", limit = c(0.25, 150)) +
    scale_y_continuous(trans = "logit", limit = c(0.875, 0.9981),
		       breaks = c(.85, .90, .95, .99, .995, .998)) +
    geom_label(size = 3, show.legend = FALSE)

falac-ecolo0.png

Figure 11: Relación lineal en los grupos

Las nuevas opciones aquí son limit, que permite variar el rango del eje, logit y breaks. Breaks permite personalizar dónde aparecen las marcas del eje, nada especialmente complicado. Logit se define como el logaritmo de la 'oportunidad relativa', 'razón de momios' u otros términos. Viene dado por la expresión log(p/)1-p)) siendo p una probabilidad o proporción. Es decir, si p es la probabilidad de supervivencia infantil, la oportunidad relativa p/(1-p) es la relación entre los que sobreviven y los que no. El aplicar un logaritmo hace que las cifras de oportunidad relativa sean simétricas, de manera que es 0 cuando las probabilidades son las mismas para cada caso.

Esto es especialmente útil cuando las escalas están muy cerca de 0 o de 1, como es el caso: en los peores casos, la tasa de mortalidad está en torno al 6% (0.94) y en los mejores en torno al 0.5% (0.995). Esta relación entre los parámetros en el numerador y el denominador hace que las variaciones se noten muchísimo (99.9 /0.1 es 10 veces superior (999) a 99/1 (99) y a su vez 10 veces superior a 90/10, que es 9) y el logaritmo hace que los incrementos sean constantes.

Todo esto es para mostrar que efectivamente hay una relación lineal entre ingresos y tasa de supervivencia infantil, pero darla por cierta fuera de estos valores medios por continente, es lo que se denomina la falacia ecológica. Si se ven los datos por país, está claro que la variabilidad es muy alta y, por ejemplo, Angola y Croacia, a similar renta per cápita (en torno a los 8$ diarios), tienen tasas de supervivencia infantil muy diferentes (alrededor de 90% en Angola y más de 99.5% en Croacia, lo que significa que un niño en Croacia tiene más o menos 20 veces más probabilidades de sobrevivir que en Angola).

falac-ecolo1.png

Figure 12: Relaciones por país

3.6 Principios de visualización de datos

3.6.1 Introducción

La mayoría de lo expuesto a continuación está basado en una charla de Karl Broman titulada «Creating Effective Figures and Tables» y en las notas de Peter Aldhous «Introduction to Data Visualization».

Los principios básicos están basados en cómo los humanos hacemos comparaciones y usamos nuestras capacidades visuales. Se compararán gráficos que siguen estos principios con otros que no lo hacen para extraer generalizaciones útiles.

3.6.2 Codificar datos usando notas visuales

Los gráficos de sectores (pie chart) muestran los datos usando áreas y ángulos proporcionales, pero los estudios sobre percepción muestran que los humanos no somos muy buenos con estos parámetros, por lo que no son formas muy recomendables de representar datos por más que los haya popularizado Excel. Los dónuts son aún peores, porque su capacidad de codificar datos solo reside en el área de las secciones del anillo. En caso de necesitar un gráfico de este estilo, se hace imprescindible acompañar la imagen de los porcentajes escritos sobre los sectores. Es mucho más recomendable una simple tabla con los valores, por ejemplo.

Los mejores gráficos para representar porcentajes son los basados en longitud y posición, como los gráficos de barras. Son fáciles de leer a simple vista y de comparar cuando se ponen uno al lado del otro.

El color y el brillo son aún peores que los ángulos y el área cuando se trata de representar datos. Aún así, pueden resultar útiles si hay que codificar más de dos dimensiones en un mismo gráfico.

3.6.3 ¿Cuándo incluir el 0?

Cuando se hace un gráfico de barras, no es correcto no comenzar las barras en 0. Si se comienza en otro valor, las diferencias entre barras parecen mucho mayores que si se hace como debe hacerse. Ejemplo: entre 16,000 y 17,000 hay un 6% de diferencia, por lo que las barras tendrán un 6% de diferencia en la longitud. Pero si el origen de las barras es 15,000, la segunda parece el doble de la primera. Es una táctica muy usada en política y prensa para presentar datos de forma tendenciosa.

Sin embargo, cuando es solo la posición y no la longitud lo que codifica los datos, como en un gráfico de dispersión, puede no incluirse el 0. De hecho, suele ser conveniente pues permite ajustar el área del gráfico al rango de datos y hacer que los puntos no se apiñen en una zona dejando espacio vacío inútil; así están más separados unos de otros y la información es más clara.

3.6.4 No distorsionar las cantidades

Suena obvio, pero es fácil (y conveniente si hay intención de manipular la información) que suceda si se escoge un mal parámetro para codificar la información, como el radio de un círculo en lugar de su área. De esta manera, el área se ve influenciada con el cuadrado del radio y las diferencias se magnifican muchísimo. Hay que limitarse al área (es el comportamiento por defecto de ggplot) o mejor, usar un gráfico de barras.

3.6.5 Ordenar los datos por valores significativos

ggplot siempre ordena las categorías en un gráfico de barras (si son alfanuméricas o factores) alfabéticamente. Pero el orden alfabético es arbitrario, y las categorías deben ordenarse de mayor a menor (o viceversa), y mantenerse este orden en el gráfico contiguo en caso de que sean comparables (como dos puntos en el tiempo, por ejemplo).

Para ordenar convenientemente los elementos de un gráfico se usa la función reorder. Más arriba hay ejemplos del uso de esta función.

3.6.6 Mostrar los datos

En muchos ejemplos anteriores se comparan datos de distintas categorías, pero cuando se intentan mostrar los resultados por grupos algunos tipos de gráficos popularizados por MS Excel, sobre todo, no son muy adecuados. Por ejemplo, los diagramas de barras que muestran la media (con o sin errores estándar añadidos encima) van desde 0 hasta la media, si lo que se representa es la altura media de un grupo ¿significa que hay humanos que mides 1cm o 2cm?

Simplemente mostrar los datos puede ser mucho más efectivo que un gráfico tan pobre en información (ver 10.5, pág. 209 en el libro para una imagen de ejemplo). El siguiente código muestra todos los datos de las alturas de hombres y mujeres del conjunto de datos heights y da una idea de cómo se distribuyen los datos, el rango para cada grupo, el valor más común, etc:

heights %>% ggplot(aes(sex, height)) + geom_point()

Aún así, como muchos de los puntos se superponen, no se puede ver toda la información. Una manera de mejorar este gráfico es añadir dos capas más: jitter 'temblor' 'fluctuación', que consiste en añadir una pequeña separación horizontal a los puntos y que no estén todos en la misma línea; y alpha blending o transparencia, que consiste en que los puntos que se superponen generan un área de mayor densidad de color, más oscura. Ver imágenes de ejemplo en el libro.

Sin embargo, el modo más efectivo de mostrar los datos es sin duda mostrar la distribución.

3.6.7 Estandarizar los ejes cuando se comparen datos

Si los ejes de dos gráficos (por ejemplo, histograma de la distribución de alturas en hombres y mujeres) no son iguales, no es fácil apreciar a simple vista que los hombres son, en promedio, ligeramente más altos. Hay que ver el rango de valores de los ejes y dónde se sitúa la media, el valor máximo, etc. Sin embargo, cuando los ejes son iguales, el desplazamiento hacia la izquierda o derecha de la forma de la distribución aporta gran información al primer golpe de vista.

Al mismo tiempo, es conveniente agrupar los gráficos verticalmente u horizontalmente dependiendo de qué tipo de dato se quiere comparar. Si los datos se diferencian horizontalmente (como en el ejemplo anterior, en el que la diferencia de altura hace que la campana se desplace hacia la derecha o la izquierda), apilar los gráficos verticalmente hace más sencillo ver las diferencias a primera vista. Si los datos se diferencian verticalmente será mejor ponerlos horizontalmente uno al lado del otro.

Ejemplos en el libro, apartado 10.6.2.

3.6.8 Considerar la posibilidad de transformaciones

Transformaciones logarítmicas, logit o la raíz cuadrada pueden ser útiles para presentar los datos de manera lógica y que se comprendan. Ver los ejemplos del apartado 10.6.3 del libro para entender cómo el uso de la media en un gráfico de barras no da una información útil y cómo la transformación logarítmica del eje Y unido a un diagrama de cajas y bigotes proporciona mucha más información y más clara.

3.6.9 Los elementos comparados deben ser adyacentes

El orden alfabético o agrupado por el componente incorrecto puede ser causa de que el gráfico no sea claro. Por ejemplo (ver 10.6.4), para comparar datos de continentes en dos años diferentes, lo óptimo es juntar las barras o diagramas que se usen de dos en dos, uno para cada año. Y hacer lo mismo para cada continente. Si, aún encima, se usan colores, todo será mucho más claro de un vistazo, que es lo que se pretende.

3.6.10 Usar colores

Usar colores lo hace todo más fácil. ggplot permite variar fácilmente los colores por defecto usando códigos RGB. Además, los colores por defecto no son buenos si se piensa en los daltónicos.

Para usar colores diferentes se define primero una paleta y luego se añade una capa con la paleta que se quiere usar al código del gráfico:

paleta_de_colores <- c("#23c1c1", "#999999", "#ff3300")
gráfico + scale_color_manual(values = paleta_de_colores)

Aquí, gráfico es el código completo del gráfico, con su filter, su ggplot, su geompoint y todo. Solo se añade una capa a mayores.

3.6.11 Gráficos de pendientes slope charts

Hasta ahora se han comparado y descrito gráficos en los que se estudia la relación entre dos variables (mortalidad e ingresos, altura y sexo, etc.), y un diagrama de densidad, scatter plot, es lo mejor y más claro. Sin embargo, hay casos en los que usar posición y ángulo para mostrar diferencias puede ser una buena idea.

Por ejemplo, cuando los puntos a mostrar no son muchos y lo que se pretende poner de manifiesto es una evolución temporal. El gráfico de la página 218 del libro se construye con el siguiente código, pues lamentablemente ggplot no tiene una manera predefinida de generar un gráfico de pendientes y hay que usar geomline para eso.

west <- c("Western Europe","Northern Europe","Southern Europe",
	    "Northern America","Australia and New Zealand")

dat <- gapminder %>%
  filter(year%in% c(2010, 2015) & region %in% west &
	   !is.na(life_expectancy) & population > 10^7)

dat %>%
  mutate(location = ifelse(year == 2010, 1, 2),
	 location = ifelse(year == 2015 &
	   country %in% c("United Kingdom", "Portugal"),
	   location+0.22, location),
	 hjust = ifelse(year == 2010, 1, 0)) %>%
  mutate(year = as.factor(year)) %>%
  ggplot(aes(year, life_expectancy, group = country)) +
    geom_line(aes(color = country), show.legend = FALSE) +
    geom_text(aes(x = location, label = country, hjust = hjust),
    show.legend = FALSE) +
    xlab("") + ylab("Life Expectancy")

Estos son los cambios en la esperanza de vida de los países occidentales con población mayor a 10 millones de habitantes.

Se definen los que son los países occidentales como los pertenecientes a las regiones especificadas en el vector west, y luego se filtran los datos de dos años específicamente y población mayor a 10 millones.

Luego se introduce una nueva variable usando mutate llamada location a la que se le da el valor 1 si el año es 2010 y 2 si no lo es (lo que implica que es 2015, pues ya hubo un filtro previo por año). Además, location será location + 0.22 si los países son UK o Portugal y el año es 2015. Esto es para desplazar las etiquetas de estos países según se ve en la imagen.

Se justifica el texto horizontalmente para el año 2010, es decir, la serie de la izquierda. Y, finalmente, se dibuja el gráfico con una capa geomline que aplica diferentes colores a los diferentes países.

La última línea lo que hace es cambiar las etiquetas (labels) de los ejes X e Y por nada y Life Expectancy, respectivamente.

Así es fácil comparar dos series temporales. Hacer lo mismo con una distribución de densidad o scatter plot es más complicado. No el código, que es más sencillo, sino la interpretación: habría que mirar el valor del eje X y compararlo con el valor del eje Y, 2010 y 2015 respectivamente. La información está ahí, pero es más difícil de extraer y no se hace al primer golpe de vista.

3.6.12 Gráfico de Bland-Altman

Es otro tipo de gráfico para estudiar diferencias, también conocido como gráfico de medias y diferencias de Tukey o gráfico MA. Ya que se muestran diferencias, se dedica un eje (el eje Y) a las diferencias entre ambas series. Y en el eje X se expresa la media de ambas medidas, lo que da una idea general. Hay ejemplo en la página 219.

3.6.13 Codificar una tercera variable

En los diagramas de dispersión vistos hasta ahora se relacionan dos variables (tasa de mortalidad e ingresos, por ejemplo). Pero se pueden codificar más variables en varias características de los puntos correspondientes a los datos, como son forma (círculos, cuadrados y triángulos), color (ya visto y usado para codificar la región o el continente) o tamaño (círculos o cuadrados más o menos grandes según el valor de la variable -ideal para variables continuas-).

Para el ajuste de las formas se usa el parámetro shape (las formas disponibles están en un cuadro en el libro, página 220, junto a un ejemplo de gráfico combinando estas codificaciones de variables).

3.6.14 Caso de estudio: Vacunas

El objetivo de este ejercicio es ver cómo representar tres variables en un solo gráfico, emulando el gráfico sobre la eficacia de las vacunas con el ejemplo del sarampión publicado por el Wall Street Journal:

533a.png

Figure 13: Gráfico original

En el gráfico se representa el tiempo en el eje X, los distintos estados en el eje Y y el ratio de casos está codificado en los colores de los cuadrados (un cuadrado por estado y año). La verdad es que es muy llamativo, pero se puede mejorar. En concreto, hay que tener en cuenta que los colores usados son buenos para casos en los que hay un centro con sentido en la escala de datos y la escala se mueve tanto hacia arriba como hacia abajo. Pero para un caso como este en el que hay una variable continua, lo mejor es una escala monocromática, en la que hay una clara conciencia de lo que es un valor más alto (color más oscuro) y un valor más bajo (color más claro)19.

Para llegar a este gráfico hay que hacer varias cosas:

  1. Como el intervalo de tiempo empieza en los años 30, hay que filtrar los estados de Alaska y Hawaii (se convirtieron en estados en los 60).
  2. Extraer los datos solo de la enfermedad que nos concierne, el sarampión (en inglés, measles).
  3. Calcular la tasa de casos por cada 10.000 habitantes y añadirla como una nueva variable en el conjunto de datos.

El código necesario para esta primera parte es el siguiente:

data(us_contagious_diseases)
la_enfermedad <- "Measles"
datos <- us_contagious_diseases %>%
  filter(!state %in% c("Alaska", "Hawaii") & disease == la_enfermedad) %>%
  mutate(rate = count / population * 10000) %>%
  mutate(state = reorder(state, rate))

En realidad, el primer paso es innecesario, podría perfectamente eliminarse la primera línea y añadir el nombre en la segunda: & disease = "Measles") %>%=. Ahora, con esta información es sencillo hacer un gráfico para ver la evolución de los casos en un estado, por ejemplo, California:

datos %>% filter (state == "California") %>%
  ggplot(aes(year, rate)) +
  geom_line() + ylab("Casos por 10.000 habitantes") +
  geom_vline(xintercept = 1963, col = "blue")

La última línea lo que hace es añadir una línea vertical (vline) que intercepta el eje X en el valor 1963 para marcar el momento de introducción de la vacuna.

Para representar tres variables en el mismo gráfico se usarán, igual que en el ejemplo del WSJ, el eje X para el tiempo, el eje Y para los estados y los colores para la tasa de aparición de la enfermedad. Se usa la geometría geom_tile para «embaldosar» el gráfico con los colores que representan la variable rate (tile 'baldosa'). También hay que incluir la transformación de raíz cuadrada para evitar que los valores extremadamente altos enmascaren la información útil del resto de valores. Recuérdese que datos ya había sido definido antes, lo único que se cambia ahora es el gráfico, de geom_line a geom_tile. A continuación el código completo:

# Cargar las bibliotecas necesarias
library(dslabs)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)

#Cargar el conjunto de datos
data(us_contagious_diseases)

# Definir la enfermedad a representar
la_enfermedad  <- "Measles"

# Filtrar los datos y guardarlos en un elemento. Añadir la tasa de incidencia de la enfermedad como una variable más y reordenar los estados.
datos <- us_contagious_diseases %>%
    filter(!state %in% c("Alaska", "Hawaii") & disease == la_enfermedad) %>%
    mutate(rate = count/population * 10000) %>%
    mutate(state = reorder(state, rate))

# Guardar el gráfico como png
png(file = "./incl/533b.png",
    width = 640,
    height = 640)

# Código del gráfico
datos %>% ggplot(aes(year, state, fill = rate)) +
    geom_tile(color = "grey50") +
    scale_x_continuous(expand = c(0,0)) +
    scale_fill_gradientn(colors = brewer.pal(9, "Reds"), trans = "sqrt") +
    geom_vline(xintercept = 1963, col = "blue") +
    theme_minimal() + theme(panel.grid = element_blank()) +
    ggtitle("Sarampión") +
    xlab("") +
    ylab("")

# Devolver el control a la consola después de guardar el gráfico
dev.off()

Después de cargar datos y bibliotecas, y filtrar los datos (igual que más arriba), la tercera línea lo que hace es cambiar el comportamiento por defecto de scale_x_continuous, que añade un 5% a partir del valor máximo en las variables continuas para acomodar bien los datos sin que lleguen al borde del gráfico. Pero en este caso, lo interesante es que las «baldosas» lleguen al borde del gráfico y lo llenen todo, de manera que se le indica que se expanda cero.

El resultado es bastante impresionante: #+caption 533b.png

Sin embargo, se usan los colores para representar valores (el número de casos de la enfermedad se representa con los distintos tonos de rojo), cosa muy poco recomendable por poco exacta. Lo mejor, como se dijo, es usar longitud y posición (barras y líneas, por ejemplo). Por eso, el siguiente ejemplo es no tan vistoso pero más informativo.

Se asigna el mismo color a los 50 estados, representados por una línea diferente. Una línea más gruesa marca la media para todo el país.

El código no me funciona, y lo pondré aquí en cuanto entienda qué pasa.

library(dslabs)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
data(us_contagious_diseases)

la_enfermedad  <- "Measles"

datos <- us_contagious_diseases %>%
    filter(!state %in% c("Alaska", "Hawaii") & disease == la_enfermedad) %>%
    mutate(rate = count/population * 10000) %>%
    mutate(state = reorder(state, rate))

		 # Calcular la media para los EEUU.

media <- us_contagious_diseases %>%
    filter(disease == la_enfermedad) %>%
    group_by(year) %>%
    summarize(tasa_en_usa = sum(count, na.rm = TRUE)/sum(population, na.rm = TRUE)*10000)

	   # Gráfico de líneas para todos los estados

png(file = "./533c.png",
    width = 640,
    height = 640)

datos %>% filter(!is.na(rate)) %>%
    ggplot() +
    geom_line(aes(year, rate, group = state), color = "grey50", 
	      show.legend = FALSE, alpha = 0.2, size = 1) +
    geom_line(mapping = aes(year, tasa_en_usa), data = media,
	      size = 1, col = "black") +
    scale_y_continuous(trans = "sqrt", breaks = c(5, 25, 125, 300)) +
    ggtitle("Casos por cada 10,000 habitantes, por estado") +
    xlab("") +
    ylab("") +
    geom_text(data = data.frame(x = 1955, y = 50),
	mapping = aes(x, y, label = "Media de los EEUU"), color = "black") +
    geom_vline(xintercept = 1963, col = "blue")

dev.off()

3.6.15 No usar gráficos 3D y tener cuidado con los dígitos significativos

Es poco informativo (nada informativo) y visualmente confuso el usar gráficos en los que se hace una representación 3D de los valores. Lo mejor es usar colores para distinguir diferentes variables en líneas o barras, pero siempre en gráficos planos (2D), que es como son representados en pantallas o papel.

En cuanto a los digitos significativos, es poco útil usar los siete ¡7! dígitos que muestra R por defecto. Hacen las tablas difíciles de leer y no aportan información importante. Lo mejor es restringir a uno o dos como mucho (salvo en casos en los que sea algo útil de verdad). Para eso se pueden usar las funciones signif o round. O, también, establecer una opción que afecte a todo el código:

options(digits = n)

donde n es el número de dígitos que se desean20.

Por último, es mejor que los datos que se desean comparar en una tabla estén dispuestos en columnas en lugar de en filas. Aunque se pueden comparar exactamente igual, se hace más difícil comparar cifras todas en la misma fila.

4 Probabilidad

4.1 Probabilidad discreta

;properties; :customid: subsec:proba-discr

4.1.1 Introducción

Las teorías acerca de la probabilidad provienen de los estudios matemáticos de los juegos de azar llevados a cabo por numerosos estudiosos, incluidos famosos matemáticos. Estas teorías se usan en los juegos de azar (obvio), en los casinos (en los que se estudian los juegos para tener beneficios) y, y esto es lo que nos interesa, en los conjuntos de datos en los que el azar tiene algún efecto.

La definición de probabilidad es muy intuitiva cuando se trata de los juegos de azar: la probabilidad de sacar un 6 al tirar un dado es de 1 entre 6. Además, el uso cotidiano de la palabra es mucho más general y amplio, menos exacto, como en la probabilidad de tener gemelos, de ser alcanzado por un rayo o la probabilidad de que X gane las elecciones.21 Estos usos de la probabilidad se basan en la inferencia estadística, cosa que se verá en el siguiente capítulo. Pero la inferencia estadística está fuertemente asentada en la probabilidad, así que hay que pasar antes por aquí.

La motivación de este capítulo son las circunstancias alrededor de la crisis financiera de 2007-2008, causada (entre otras cosas) por subestimar ciertos riesgos vendidos por las instituciones financieras. Las hipotecas y otros productos fueron vendidos a un precio que asumía que los trabajadores pagarían sus deudas mes a mes, pero esto no resultó así por una combinación de factores y los precios se hundieron. El resto es historia.

Los conceptos más importantes que se van a introducir son:

  • Variables aleatorias.
  • Independencia.
  • Simulaciones Monte Carlo.
  • Error estándar.
  • Margen de error.
  • Valores esperados.
  • Teorema del valor central.

4.1.2 Introducción a la probabilidad discreta

Corresponde al capítulo 13 del libro de texto.

La probabilidad discreta se refiere a las variables categóricas, aquellas en las que los elementos pueden tomar un número reducido de valores (las caras de un dado, los estados de EEUU, las cartas de la baraja, etc.). Más tarde se introducirán las variables continuas, más habituales en la ciencia de datos.

Como se comentó en la introducción, la palabra probabilidad tiene un uso coloquial muy habitual, como ¿Cuál es la probabilidad de que llueva mañana? Es muy difícil responder preguntas acerca de la probabiliad sin una definición más exacta y matemática de lo que es la probabilidad y a qué se ciñe. Por ejemplo, en casos como si hay tres canicas azules y dos rojas en una urna y saco una al azar, cuál es la probabilidad de que sea roja? Intuitivamente es un 40% o 2/5 (dos entre cinco), pero hace falta una definición precisa para poder extrapolar el ejemplo a casos más complejos.

En este caso de ejemplo (las canicas en la urna), hay dos resultados de cinco posibles que satisfacen la condición necesaria para el evento «sacar una canica roja». Es decir, dos entre cinco o 2/5 = 0.4 = 40%. Y esto pasa porque todos los resultados tienen la misma probabilidad de suceder cuando se saca una canica de la urna.

La probabilidad de un suceso es la proporción de veces que el suceso acontece si se repite el experimento una y otra vez, independientemente y en las mismas condiciones.

Y un evento es

algo que puede suceder cuando el caso que se trata sucede por azar.

En el ejemplo anterior, el suceso o evento es «sacar una bola roja». En otros casos, puede ser sacar 42 votantes del partido A y 58 del partido B de entre 100 votantes escogidos al azar. En variables continuas, un suceso o evento será algo diferente, como «la persona es más alta de 1.60m». Este uso de la probabilidad se vio antes en el apartado referido a las distribuciones normales.

La notación de aquí en adelante será:

Pr(A) : probabilidad de que suceda el evento A. x > 1.60: para variables continuas.

  1. Simulaciones Monte Carlo

    En R hay funciones que permiten simular el acto de escoger elementos al azar dentro de un conjunto dado, como en el ejemplo anterior de sacar una bola, roja o azul, de entre un conjunto de 5 bolas metidas en una urna. Una de estas es la función sample.

    Para hacer una simulación de una urna de la que se extraerán bolas al azar, primero hay que generar la urna, es decir, generar un conjunto en el que haya 5 bolas, 3 azules y 2 rojas. La función rep() replica un elemento un cierto número de veces. El siguiente código es equivalente a escribir los colores de las bolas uno por uno22:

    bolas <- rep(c("azul", "roja"), times = c(3, 2))
    

    Dado que times es un argumento obligado y que aparece en segundo lugar, escribir rep(c("azul", "roja"), c(3, 2)) es exactamente lo mismo.

    Ya está configurada la urna en la que están contenidas las cinco bolas, ahora hay que extraer una bola de ahi:

    sample(bolas, 1)
    

    Es decir, extraer un elemento al azar del conjunto bolas. Para lograr una estadística fiable, este experimento de extraer una bola del conjunto debería ser repetido continuamente, cosa que no es posible hacer, pero sí un gran número de veces, lo que permite hacernos una idea de cómo se aproximan los números a lo que la teoría dice. Más tarde se verá cómo decidir lo que es «un número suficiente de veces». La función replicate lo que hace es repetir una operación un cierto número de veces.

    B <- 10000
    sucesos <- replicate(B, sample(bolas, 1))
    tabla <- table(sucesos)
    tabla
    

    Explicación del código: se asocia un número a un elemento, solo para no escribirlo directamente en la función. Luego se asocia la función de replicar la extracción de bolas a un elemento sucesos, del cual posteriormente se crea una tabla, la cual se muestra. Una opción es ver las proporciones en lugar de la cuenta del número de resultados: prop.table(tabla).

    Aunque los números no sean exactamente 0.6 y 0.4, son muy aproximados, y haciendo B más grande, los valores se irán aproximando paulatinamente a estos límites. Estos están definidos cuando B es infinito, es decir, el experimento se repite para siempre. Además, son muy fáciles de calcular teóricamente, incluso intuitivamente, como se vio más arriba. Sin embargo, hay casos en que no es fácil (o no es posible) tener una visión teórica de cuál debe ser la probabilidad, y las simulaciones Monte Carlo proporcionan una vista de los valores aproximados.

    La función sample extrae el número de elementos que se le pasa como parámetro sin reemplazo, o sea, en la urna quedan cuatro bolas tras extraer una, quedan tres tras extraer dos, etc. De hecho, si se pide una extracción de seis bolas sample(bolas, 6), devolverá un error porque no hay elementos suficientes. Pero hay un argumento, replace, que está desactivado por defecto, que puede activar el reemplazo, esto es, devolver la bola a la urna tras extrerla y comprobar el color. De este modo, se puede hacer la simulación Monte Carlo sin la función replicate, simplemente sacando B bolas de la urna reemplazándolas cada vez.

    sucesos <- sample(bolas, B, replace = TRUE)
    prop.table(table(sucesos))
    

    Ver una explicación más detallada de cómo funciona una simulación Monte Carlo más abajo.

    Otro ejemplo de simulación Monte Carlo (usando el ejemplo de más abajo para explicar la independencia):

    El equipo A juega una liga de cuatro partidos contra el equipo B. El equipo B es mejor y tiene un 60% de probabilidades de ganar los partidos que dispute contra A (cada uno de ellos). Por tanto, dado que son sucesos independientes, la probabilidad de que B gane los cuatro partidos es de 0.6^4. ¿Cuál es la probabilidad de que A gane al menos uno de los partidos? La respuesta es que, como son suposiciones opuestas, 1-(0.6^4) o 0.87.

    Para escribir un código que simule un número grande de ligas de cuatro partidos y ver cuál es el resultado:

        # La variable B especifica el número de veces que se
        # repite la simulación.
    B <- 10000
        # 'gana_a' es un vector lógico creado por la repetición de la
        # operación de 'liga_simulada'. El resultado de comprobar si hay algun
        # 'gana' en la operación se almacena en el vector lógico, y así B
        # veces.
    gana_a  <- replicate(B, {
        liga_simulada <- sample(c("pierde","gana"), 4, replace = TRUE,
    			       prob = c(0.6, 0.4))
        any(liga_simulada %in% "gana")
       })
        # Calcular la frecuencia de B iteraciones en las que el equipo A gana
        # al menos un partido -cuántos TRUE hay en el vector.
    mean(gana_a)
    
  2. Configurar el valor aleatorio semilla

    En los ejemplos y en el libro de texto se usan generadores de números aleatorios, lo cual hace que cada vez que se repite el experimento los resultados pueden diferir. No pasa nada porque ese es precisamente el comportamiento esperado de los números aleatorios.

    Sin embargo, por cuestiones de reproducibilidad, es posible «ajustar» la semilla del generador de números aleatorios (RNG, Random Numbers Generator). En el libro y los ejemplos del curso está ajustada al valor 1986.

    set.seed(1986)
    

    Es un modo habitual de escoger este valor, restar al año actual el mes actual y el día actual, de manera que 1986 = 2018-12-20 porque se ajustó el 20 de diciembre de 2018.

    Además, con la actualización a R 3.6 cambió el modo de ajustar este «valor semilla». Para reproducir el comportamiento original, son necesarios ciertos parámetros:

    set.seed(1, sample.kind = "Rounding")
    
  3. Usar la función mean para calcular probabilidades

    La función mean, que calcula la media de una serie de valores, cuando se aplica a un vector lógico, muestra la proporción de TRUE en ese vector. Si se aplica al anterior vector bolas, se puede ver la proporción de cada uno de los factores, es decir, la probabilidad de que salga escogido cada uno de ellos:

    mean(bolas == "azul")
    [1] 0.6
    

    El funcionamiento de este proceso consta de varios pasos: primero se genera un vector lógico intermedio de elementos TRUE y elementos FALSE dependiendo de si cada uno de los elementos del vector bolas es igual a azul o no. Luego, al aplicar la función mean, el vector intermedio es forzado a un vector numérico en el que los TRUE son convertidos en unos (1) y los FALSE son convertidos en ceros (0). Finalmente, se calcula la media (número de unos dividido entre el número de elementos, tanto unos como ceros).

    En datos discretos, como las bolas en una urna o la proporción de votantes de un cierto número (limitado) de partidos políticos, las cartas de la baraja o las caras de un dado, las probabilidades de que suceda un evento (se extraiga una bola de un color determinado, el dado muestre un número concreto, etc.) vienen dadas por la proporción de elementos de esa clase contenidos en el conjunto.

    En variables categóricas discretas, la distribución de probabilidades viene dada por las proporciones de elementos de cada grupo en el conjunto de datos.

  4. Independencia

    Dos sucesos son independientes si el resultado de uno no afecta al resultado del otro.

    Un ejemplo clásico son las dos caras de una moneda: no importa el resultado cuando se tira una moneda y sale cara o cruz, al volver a tirarla la probabilidad sigue siendo el 50%. Otro ejemplo es sacar las bolas de la urna con reemplazo. Al sacar la segunda bola, las probabilidades están intactas sin importar el color de la primera extracción.

    Sucesos no independientes son los que suceden con las barajas de cartas. Las probabilidades de sacar un rey de una pila de cartas es de 4 entre 52, o lo que es lo mismo, 1 entre 13 (1/13 = 7.7%)23. Es lo mismo que antes. Pero si la carta no vuelve al mazo, las probabilidades de sacar un rey no son 4/52 sino 3/51, pues hay una carta menos y un rey menos en el mazo. Si la primera carta no es un rey, las probabilidades son 4/51. Por lo tanto, ambos sucesos no son independientes, el resultado de uno afecta al resultado del otro.

    Un caso extremo de sucesos no independientes es sacar bolas de la urna sin reemplazarlas. Si se extraen cinco bolas y no se conocen los resultados de las extracciones (por ejemplo, asignándolas a un vector y no comprobando el contenido del vector: x <- sample(bolas, 5)), se puede intuir que el resultado de la primera extracción es de 40% de probabilidades de que la bola sea roja y 60% de que sea azul. Pero si se comprueban las extracciones 2 a 5: x[2:5] y los resultados muestran dos bolas rojas, las suposiciones sobre la primera extracción cambian completamente; ahora hay un 100% de probabilidades de que la primera bola sea azul. Si fuesen tiradas de una moneda o las bolas fueran devueltas a la urna, no sería igual.

    Esto se denomina probabilidad condicional: la probabilidad de que un suceso B ocurra dado el acontecimiento del suceso A. Los cálculos condicionales son necesarios para llegar a conclusiones correctas. La notación que se utiliza es la siguiente:

    Pr(Carta 2 es rey | Carta 1 es rey) = 3/51

    Es decir, la probabilidad de que la segunda carta sea rey con la condición de que la primera carta sea rey es 3/51. En general, los sucesos independientes siguen la siguiente fórmula:

    Pr(A | B) = Pr(A)

    O lo que es lo mismo, la probabilidad de que suceda A (sacar cruz) dado que antes ha sucedido B (sacar cara) es la misma que la de que suceda A sin tener en cuenta lo que haya sucedido antes.

    Para conocer la probabilidad de que dos sucesos A y B sucedan, se usa la regla de la multiplicación (para dos y tres sucesos dependientes):

    Pr(A y B) = Pr(A) * Pr(B | A) Pr(A y B y C) = Pr(A) * Pr(B | A) * Pr(C | A y B)

    Y para tres sucesos independientes (mucho más sencillo):

    Pr (A y B y C) = Pr(A) * Pr(B) * Pr(C)

    Ahora bien, es sumamente complicado identificar la dependencia o independencia de dos sucesos fuera de las cartas, las bolas en la urna o las caras de la moneda. Un ejemplo práctico (extraído de la realidad): en un tribunal se juzga un caso en el que se sabe que el sospechoso es un hombre con bigote y barba. Si el 10% de los hombres tiene barba y, al mismo tiempo, el 20% de los hombres tiene bigote, la regla de la multiplicación dice que hay un 2% de probabilidades de que un hombre tenga barba y bigote. Pero tal asunción no es cierta, el hecho de tener barba influye mucho en el hecho de tener bigote, hasta el punto de que Pr(bigote | barba) = 0.95. Por lo tanto, el cálculo correcto sería 1/10 * 9.5/10 = 9.5%.

    Para hacer los cálculos usando un código en R, se puede hacer lo siguiente. Supongamos que en la urna hay tres bolas azules, 5 bolas rojas y 7 bolas amarillas (el total es 15). La probabilidad de que la primera extracción sea azul es

    azul <- 3
    roja <- 5
    amarilla <- 7
    
    p_1 <- azul / (azul + roja + amarilla)
    
    p_1
    

    Así de simple. La probabilidad de que la segunda extracción no sea azul dependerá de si hay o no hay reemplazo, es decir, de si los sucesos son dependientes o independientes. Si no hay reemplazo

    azul <- 3
    roja <- 5
    amarilla <- 7
    
    p_1 <- azul / (azul + roja + amarilla)
    p_2 <- (roja + amarilla) / ((azul - 1) + roja + amarilla)
    
    p_1 * p_2
    

    Y si hay reemplazo (es decir, ambos sucesos son independientes):

    azul <- 3
    roja <- 5
    amarilla <- 7
    
    p_1 <- azul / (azul + roja + amarilla)
    p_2 <- (roja + amarilla) / (azul + roja + amarilla)
    
    p_1 * p_2
    

    Un ejemplo con un dado. Las tiradas de dados son sucesos independientes porque no se puede «eliminar» una cara de un dado una vez que ha salido, cada vez que se tira hay seis posibles resultados con la misma probabilidad de salir.

    La probabilidad de que no salga un 6 cuando se tira un dado es de 5/6 (de 6 posibles resultados, 5 cumplen la condición de no ser un 6). Sucesivas tiradas son sucesos completamente independientes de este. Por lo tanto, ¿cuál es la probabilidad de que no salga un 6 en seis tiradas sucesivas? Según la regla de la multiplicación de más arriba, se multiplica la probabilidad de cada uno de los sucesos independientes, que en este caso es la misma todas las veces: 5/6 * 5/6 * 5/6 * 5/6 * 5/6 * 5/6 o, lo que es lo mismo, 5/6^6.

    prob_6 <- 5/6
    prob_6_veces <- prob_6 ^ 6
    

    Hay que recordar que, cuando se calculan probabilidades, puede ser interesante calcular la probabilidad contraria y restarla de 1. Por eso, es mejor recordar que lo contrario de Todo es No todo o Al menos uno no. A su vez, lo contrario de Nada es Algo o Al menos uno. Por ejemplo:

    El equipo A juega una liga de cuatro partidos contra el equipo B. El equipo B es mejor y tiene un 60% de probabilidades de ganar los partidos que dispute contra A (cada uno de ellos). Por tanto, dado que son sucesos independientes, la probabilidad de que B gane los cuatro partidos es de 0.6^4. ¿Cuál es la probabilidad de que A gane al menos uno de los partidos? La respuesta es que, como son suposiciones opuestas, 1-(0.6^4) o 0.87.

4.1.3 Combinaciones y permutaciones

Corresponde al apartado 13.6 del libro de texto con el mismo título.

Hasta ahora se ha usado el ejemplo de la urna y las bolas de colores, pero los cálculos se pueden complicar bastante en casos más elaborados. Las proporciones en el ejemplo de la urna son fáciles de contar intuitivamente: si hay 5 posibilidades y de ellas, 3 son azules, pues la probabilidad de sacar azul es 3/5.

Para hacer cálculos más complejos es útil la baraja de cartas (como se dijo más arriba, será la baraja francesa). Y, por supuesto, lo primero es construir la baraja. Para eso hacen falta las funciones paste y expand.grid.

La función paste une dos cadenas de texto con un espacio en el medio. Por supuesto, también funciona con vectores alfanuméricos, uniendo ambos elemento a elemento.

numero <- "Tres"
palo <- "Picas"
paste(numero, palo)

dará como resultado Tres Picas.

Los elementos pueden ser vectores: paste(letters[1:5], as.character(1:5)) devuelve el resultado: "a 1" "b 2" "c 3" "d 4" "e 5".

La función expand.grid crea una tabla con las combinaciones posibles de los vectores que se le pasen como argumento. Por ejemplo, las posibles combinaciones de ropa de colores puede ser:

pant <- c("Azul", "Negro")
cams <- c("Blanca", "Azul", "Estampada")
expand.grid(Pantalones = pant, Camisas = cams)

El resultado24 será:

  Pantalones   Camisas
1       Azul    Blanca
2      Negro    Blanca
3       Azul      Azul
4      Negro      Azul
5       Azul Estampada
6      Negro Estampada

Con estos elementos ya se puede crear una baraja completa (francesa, repito):

palos <- c("Diamantes", "Tréboles", "Corazones", "Picas")
numeros <- c("As", "Dos", "Tres", "Cuatro", "Cinco",
	     "Seis", "Siete", "Ocho", "Nueve", "Diez",
	     "Sota", "Reina", "Rey")
baraja <- expand.grid(numero = numeros, palo = palos)
baraja <- paste(baraja$numero, baraja$palo)

La baraja está lista. Ahora hay que comprobar las probabilidades de ciertos sucesos, como que la probabilidad de sacar un rey a la primera es 1/13.

Para ello, se calculan las proporciones de los posibles resultados que satisfacen la condición propuesta. Se crea un vector que contenga las cuatro maneras de sacar un rey (hay un rey por palo y hay cuatro palos, luego cuatro posibilidades diferentes).

reyes <- paste("Rey", palos)
mean(baraja %in% reyes)
# [1] 0.07692308

El código une la palabra Rey con el contenido del vector palos definido con anterioridad: "Rey Diamantes" "Rey Tréboles", etc. Luego se calcula la media del vector lógico que resulta de comprobar si los elementos de baraja están en reyes (ver función %in%). El resultado es exactamente 1/13 (7.69%).

El siguiente cálculo será comprobar la probabilidad de sacar un segundo rey después de sacar un primer rey. Anteriormente se dedujo que tras sacar un rey quedan 3 reyes y 51 cartas, así que la probabilidad es 3/51. Pero para confirmarlo usando código en R, hay que comprobar las posibilidades usando combinaciones y permutaciones.

Las funciones combinations() y permutations() están disponibles en el paquete gtools, de modo que hay que instalarlo y cargar la biblioteca antes de usar estas funciones: install.packages("gtools") y library(gtools).

La función permutations(n, r) calcula, para cualquier lista de n elementos todas las posibles maneras de seleccionar r elementos. En esta lista de permutaciones el orden importa, por lo que, por ejemplo, la combinación 1 3 es diferente de 3 1 y aparecen las dos. Tampoco aparecen los pares 1 1, 2 2, etc., porque una vez que «se saca» un elemento para comprobarlo, este queda «fuera de la urna». Un ejemplo de funcionamiento de permutations() sería el siguiente código, que crea una lista de 5 teléfonos al azar de entre todos los posibles números de 7 cifras25.

install.packages("gtools")
library(gtools)
   # crear la liśta de números
todos_los_numeros <- permutations(10, 7, v = 0:9)
   # guardar el número de filas (nrs de tel) en el elemento
n <- nrow(todos_los_numeros)
   #crear un índice de las 5 filas al azar
indice <- sample(n, 5)
   # mostrar los números correspondientes al índice
   # (en cuanto al número de fila)
todos_los_numeros[indice,]

Lo interesante de este trozo de código es que permutations() permite definir la longitud del vector pero también definir su contenido. Por defecto, el contenido de un vector de longitud n es 1:n, pero en este caso se define su longitud como 10 y su contenido como 0:9. El resto se explica bien con los comentarios intercalados.

Los usos de permutations() son interesantes en cuestión de probabilidades. Por ejemplo, de cuántas maneras diferentes se pueden distribuir tres medallas (oro, plata, bronce) entre los 8 atletas de una carrera olímpica de 200m lisos: permutations(8, 3) son las permutaciones posibles de 8 elementos (los corredores) tomados de 3 en 3 (las medallas disponibles). Ese código arroja una tabla de 3x336, luego son 336 maneras diferentes de distribuir las 3 medallas entre los 8 participantes.

#install.packages("gtools")
library(gtools)
palos <- c("Diamantes", "Tréboles", "Corazones", "Picas")
numeros <- c("As", "Dos", "Tres", "Cuatro", "Cinco",
	     "Seis", "Siete", "Ocho", "Nueve", "Diez",
	     "Sota", "Reina", "Rey")
baraja <- expand.grid(numero = numeros, palo = palos)
baraja <- paste(baraja$numero, baraja$palo)

    #Se crea un vector 'reyes'
reyes <- paste("Rey", palos)
    # Probabilidad de que salga un rey (media de 'TRUES'
    #en el vector lógico compuesto por los elementos de 'baraja'
    # contenidos en 'reyes')
mean(baraja %in% reyes)
    # Calcular las posibles combinaciones de las 52 cartas puestas en
    #parejas (porque se van a sacar dos).
mano <- permutations(52, 2, v = baraja)
    # Se asocia el contenido de las combinaciones a cada una de las cartas
carta_uno <- mano[,1]
carta_dos <- mano[,2]
    # Probabilidad de que la primera carta sea un rey:
    # a) posibles combinaciones con rey en la primera carta (204)
sum(carta_uno %in% reyes)
    # b) posibles combinaciones con rey en la segunda carta divididas
    # por todas las posibles (aquí 'sum' podría ser sustituido por
    # 'mean' sin problema).
sum(carta_uno %in% reyes & carta_dos %in% reyes)/
    sum(carta_uno %in% reyes)

La última línea es una «version en R» de la regla de la multiplicación (ver el apartado Independencia): la probabilidad de B dado A es la probabilidad de A y B dividida por la probabilidad de A.

Pr(B | A) = Pr(A y B) / Pr(A)

Hasta aquí, en el ejemplo de los reyes, el orden importaba (ESTO HAY QUE MIRARLO PORQUE NO LO ENTIENDO MUY BIEN), pero si el orden no importa, como en Blackjack26, se usa la función combinations() en lugar de permutations(). La función combinations() no muestra todas las combinaciones, sino que obvia las que se pueden considerar repetidas, es decir, si muestra la combinación 2 3, no muestra 3 2. Así que para calcular la probabilidad de un Blackjack directo (as + figura), el código necesario sería:

#install.packages("gtools")
library(gtools)
					# Crear la baraja
palos <- c("Diamantes", "Tréboles", "Corazones", "Picas")
numeros <- c("As", "Dos", "Tres", "Cuatro", "Cinco", "Seis", "Siete", "Ocho", "Nueve", "Diez", "Sota", "Reina", "Rey")
baraja <- expand.grid(numero = numeros, palo = palos)
baraja <- paste(baraja$numero, baraja$palo)
					# Se crea un vector 'reyes'
reyes <- paste("Rey", palos)
					# Se crea un vector con las combinaciones con los ases
ases <- paste("As", palos)
					# Se crea un vector con las figuras, luego se crea una tabla con las figuras de cada palo, y se unen los nombres de las cartas con sus palos respectivos (igual que se hizo con el elemento 'baraja'.
figuras <- c("Diez", "Sota", "Reina", "Rey")
figuras <- expand.grid(numero = figuras, palo = palos)
figuras <- paste(figuras$numero, figuras$palo)

					# Como antes, se crean combinaciones de los 52 elementos (cartas) de dos en dos, tomando como origen de datos el vector 'baraja' y se cuentan las apariciones de figuras y ases, o ases y figuras.
mano <- combinations(52, 2, v = baraja)
mean((mano[,1] %in% ases & mano[,2] %in% figuras) |
     (mano[,1] %in% figuras & mano[,2] %in% ases))

También podría haberse hecho, en lugar de calcular las probabilidades «matemáticamente», hacer una simulación Monte Carlo para estimarlas. En este caso, una «tirada» en la que se sacan dos cartas sería mano <- sample(baraja, 2). Con esto en mente, lo que hay que hacer será repetir una y otra vez la tirada y contar los Blackjack que aparecen.

#install.packages("gtools")
library(gtools)
					# Crear la baraja
palos <- c("Diamantes", "Tréboles", "Corazones", "Picas")
numeros <- c("As", "Dos", "Tres", "Cuatro", "Cinco", "Seis", "Siete", "Ocho", "Nueve", "Diez", "Sota", "Reina", "Rey")
baraja <- expand.grid(numero = numeros, palo = palos)
baraja <- paste(baraja$numero, baraja$palo)
					# Se crea un vector 'reyes'
reyes <- paste("Rey", palos)
					# Se crea un vector con las combinaciones con los ases
ases <- paste("As", palos)
					# Se crea un vector con las figuras, luego se crea una tabla con las figuras de cada palo, y se unen los nombres de las cartas con sus palos respectivos (igual que se hizo con el elemento 'baraja'.
figuras <- c("Diez", "Sota", "Reina", "Rey")
figuras <- expand.grid(numero = figuras, palo = palos)
figuras <- paste(figuras$numero, figuras$palo)

					# Hacer una «tirada» de dos cartas y repertirla un número B de veces. Contar las veces que aparece la combinación as+figura, figura+as.
B <- 10000
resultado <- replicate(B, {
    mano <- sample(baraja, 2)
    (mano[1] %in% ases & mano[2] %in% figuras) |
	(mano[1] %in% figuras & mano[2] %in% ases)
})
mean(resultado)

A medida que B se hace mayor, aumenta la precisión de la estimación.

4.1.4 El problema del cumpleaños

Supongamos que hay un aula con 50 alumnos y queremos saber la probabilidad de que dos de ellos cumplan años el mismo día.

Un cumpleaños puede expresarse como un número entre 1 y 365, por lo que un código muy simple establece los posibles cumpleaños: 1:365. De entre estos, se extraen 50 números al azar y se comprueba si existen duplicaciones en el vector creado.

n <- 50
cumples <- sample(1:365, n, replace = TRUE)
any(duplicated(cumples))

Este código devuelve TRUE si existe algún duplicado, es decir, si hay algún número repetido en el vector cumples, lo que quiere decir que hay más de uno de los 50 alumnos con un mismo cumpleaños. Ahora solo hay que repetir la operación un gran número de veces para hacer una simulación Monte Carlo:

n <- 50
B <- 100000
    # Repeticiones del cálculo creando un vector
    # 'resultado' en el que se almacenan los TRUE
    # y FALSE de la comprobación 'any + duplicated'.
resultado <- replicate(B, {
    cumples <- sample(1:365, n, replace = TRUE)
    any(duplicated(cumples))
})
    # Comprobación del resultado
mean(resultado)

La probabilidad es realmente alta (en torno al 95%).

Sin embargo, hacer una comprobación de cuándo la probabilidad aumenta o disminuye, en qué punto la probabilidad es un valor determinado (p. ej., 50%), etc., es algo más complejo. Primero, usaremos el mismo código para crear una función que haga una simulación Monte Carlo para varios valores de n, siendo n el número de personas en el grupo. Y aquí está el problema: aunque las funciones y operaciones sobre vectores suelen aplicarse elemento a elemento, es algo que no pasa siempre y desde luego no pasa con las funciones definidas por el usuario.

    # Se crea una función 'calcula_prob' que toma como
    #parámetro obligatorio n y como parámetro por defecto
    # B = 10,000. Hace una simul. Monte Carlo y comprueba
    # si hay duplicados en el resultado

calcula_prob <- function(n, B = 10000) {
	mismo_dia <- replicate(B, {
	cumples <- sample(1:365, n, replace = TRUE)
	any(duplicated(cumples))
    })
    # La media de los TRUE de que consta el vector 'mismo_dia'
    # (resultado de la función 'any') es la probabilidad.
    mean(mismo_dia)
}
    # Se define n como una secuencia de valores de 1 a 60
    # para explorar las probabilidades con estos valores
    # del número de personas en el grupo.
n <- seq(1, 60)

Pero esto no funciona porque R no aplica calcula_prob() a todos los valores de n. Para eso hay que forzar que calcula_prob() se calcule elemento a elemento dentro del objeto n. Eso se consigue con la función sapply().27 En efecto, calcula_prob(n) devolverá 0 porque la función espera un número, no un vector. sapply() permite la aplicación de cualquier función sobre un vector, elemento a elemento.

Por tanto, hay que usar esta función pasando como parámetros el origen de datos (n en este caso) y la función a aplicar (calcula_prob). Por eso se añade el código siguiente:

probab <- sapply(n, calcula_prob)
plot(n, probab)

El resultado de aplicar calcula_prob sobre cada elemento de n se almacena en probab. Luego se puede comprobar el resultado haciendo una simple gráfica de ambos elementos.

3-123.png

Figure 14: Gráfica de la probabilidad de cumpleaños compartido

Ahora bien, no solo sería más exacto calcular la probabilidad usando la teoría y las matemáticas sino también mucho más rápido, pues no hay que repetir 10 000 veces un experimento para cada valor de n. Para simplificar al máximo, se calcula no la probabilidad de que suceda el evento sino la probabilidad de que no suceda. De modo que hay que ver la probabilidad de que los cumpleaños sean únicos, es decir, que la fecha no sea compartida.

La primera persona del grupo tiene una probabilidad del 100% de que su fecha sea única (solo es una persona, todo el año es suyo). La probabilidad de que la segunda persona tenga una fecha única es de 364/365, pues la fecha del primero ya «está cogida». La tercera persona, ya que hay dos fechas ya escogidas, tiene una probabilidad de que su fecha de cumpleaños sea única de 363/365. Y así continuamente hasta la persona 365, que tiene una probabilidad de 0 (ya no hay más días libres).

La regla de la multiplicación (ver arriba) establece que la probabilidad de que se den los sucesos A y B | A es la probabilidad de A multiplicada por la de B | A. Por eso, la probabilidad de que los n cumpleaños sean únicos es 1 · 364/365 · 363/365 · … · (365-n+1)/365. Y eso se calcula con un simple código:

    # Simulación monte Carlo
calcula_prob <- function(n, B = 10000) {
	mismo_dia <- replicate(B, {
	cumples <- sample(1:365, n, replace = TRUE)
	any(duplicated(cumples))
    })
    # La media de los TRUE de que consta el vector
    # 'mismo_dia' (resultado de la función 'any') es
    # la probabilidad.
    mean(mismo_dia)
}
    # Se define n como una secuencia de valores de 1
    # a 60 para explorar las probabilidades con estos valores
    # del número de personas en el grupo.
n <- seq(1, 60)
    # Probabilidad exacta
prob_exacta <- function(n){
    prob_unica <- seq(365, 365-n+1)/365   # vector de fracciones según fórmula
    1 - prod(prob_unica)                  # producto del vector y restarlo de 1
}
    # aplicar la función a cada n
eprob <- sapply(n, prob_exacta)

png(file = "./3-123b.png")
plot(n, prob)                   # gráfica de la simul. Monte Carlo
lines(n, eprob, col = "red")    # añadir una línea roja para la prob. exacta
dev.off()

3-123b.png

Figure 15: Monte Carlo y probabilidad exacta superpuesta

De este resultado se pueden extraer dos conclusiones importantes:

  1. La probabilidad exacta y matemática es preferible por rapidez y exactitud.
  2. Cuando no es posible, por cualquier razón, calcular la probabilidad matemáticamente, la simulación Monte Carlo es una buena aproximación con un número alto de repeticiones del experimento.

4.1.5 ¿Cuál es un buen valor de repetición para una simulación Monte Carlo?

Las pruebas hechas con las simulaciones Monte Carlo han dado muy buenos resultados con un número de repeticiones de 10 000. Sin embargo, es posible que en algunos casos este número no sea suficiente, o que sea inviable computacionalmente repetir una operación 10 000 veces para obtener un resultado y haya que escoger un número B mucho menor.

No se puede saber si la aproximación que proporciona la repetición de experimentos es fiable o no, porque no se sabe si el número es suficiente o excesivo. Hay matemáticas muy complejas detrás de la respuesta a esta pregunta, pero la aproximación práctica a continuación, en la que se observa la estabilidad de la aproximación, es muy útil.

B <- 10^seq(1, 5, len = 100)
    # Mismo código que más arriba para definir la función,
    # pero ahora se varía el valor de B y se proporciona
    # un valor fijo a n.
calcula_prob <- function(B, n = 22) {
	mismo_dia <- replicate(B, {
	cumples <- sample(1:365, n, replace = TRUE)
	any(duplicated(cumples))
    })
    mean(mismo_dia)
}
    # Calcular las probabilidades con varios valores de B.
probab <- sapply(log10(B), calcula_prob)
    # Dibujar un gráfico
png(file = "./3-124a.png")
plot(B, probab, type = "l")
dev.off()

Cuyo resultado es:

3-124a.png

Figure 16: Variación del resultado a medida que B aumenta.

4.1.6 Regla de la adición y Monty Hall

Este apartado corresponde al punto 13.5.3 del libro de texto.

  1. La regla de la adición

    En el apartado dedicado a la independencia de los eventos se explicó la regla de la multiplicación para calcular la probabilidad de dos o más sucesos, sean estos dependientes o independientes. En aquel caso, lo que se medía era la probabilidad de que se den A y B (operador AND). Ahora se medirá la probabilidad de que suceda A o B (operador OR).

    Pr (A o B) = Pr(A) + Pr(B) - Pr(A y B)

    Es decir, la probabilidad de A más la probabilidad del suceso B menos la probabilidad de que se den a la vez A y B.

    En el ejemplo del libro de texto, un blackjack directo (as y figura o figura y as), la probabilidad de que se den los dos sucesos (as y figura por un lado, figura y as por el otro) es 0, pues ambas cosas no pueden suceder (solo se sacan dos cartas). Por lo tanto, en este ejemplo concreto solo se suman ambas probabilidades. Esto es, la probabilidad de un as y luego una figura (regla de la multiplicación): 4/52 * 16/51 más la probabilidad de una figura y luego un as: 16/52 * 4/51. El resultado es el mismo que con la estimación Monte Carlo: 0.05.

  2. El problema de Monty Hall

    Este es un problema matemático basado en un concurso de televisión en el que se presentan tres puertas al concusante. Detrás de una de ellas hay un premio gordo y detrás de las otras dos un premio cutre. El concursante escoge una puerta que puede o no tener el premio gordo (digamos que es la puerta 1). El presentador28 abría entonces una de las otras dos puertas, la que no tiene el premio gordo (esta será la puerta 2). En este punto tenemos: una puerta abierta con un premio cutre y dos puertas cerradas, una escogida por el concursante y la que será la puerta 3.

    Entonces, el presentador ofrece al concursante la posibilidad de cambiar su elección. La primera impresión dice que entre la puerta 1 y la 3 (la 2 ya está descartada) la probabilidad es del 50%. ¿Cuál es la probabilidad real de escoger bien si se cambia la decisión o se mantiene?

    Se pueden estudiar las probabilidades haciendo experimentos Monte Carlo jugando 10 000 concursos con cada una de las dos tácticas.

    Código para el experimento manteniendo la puerta:

        # establecer el número de repeticiones del experimento
    
    B <- 10000
    plantarse <- replicate(B, {
        # numerar las puertas de 1 a 3 (son nombres, no números)
    	puertas <- as.character(1:3)
        # establecer los premios aleatoriamente. ¿Es una cabra un premio lo
        #  suficientemente cutre? No si eres cabrero...
    	premios <- sample(c("coche","cabra","cabra"))
        # comprobar qué puerta tiene premio: el número de puerta que
        # corresponde al elemento del vector de premios que es "coche". Es
        # decir, si "coche" es el segundo elemento del vector "premios",
        # entonces puertas[2] es la puerta premiada. Aquí se inserta una
        # condición de otro vector en lugar de un número de elemento.
    	puerta_premiada <- puertas[premios == "coche"]
        # escoger una puerta al azar
    	eleccion  <- sample(puertas, 1)
        # abrir una de las puertas sin premio que no ha sido escogida.
        # Significa: escoger 1 al azar (sample(x, 1)) dentro del vector
        # "puertas" que no sea (!puertas) ninguno de los elementos
        # correspondientes a "eleccion" ni "puerta_premiada".
    	abrir <- sample(puertas[!puertas %in% c(eleccion, puerta_premiada)],1)
        # plantarse es quedarse con la misma opción
    	plantarse <- eleccion
        # comprobar si el resultado corresponde a la puerta premiada.
    	plantarse == puerta_premiada
    })
        # probabilidad de ganar plantándose
    mean(plantarse)
    

    En realidad, la línea 12 (abrir) no es necesaria en absoluto, pues abrir o no una puerta no cambia nada ni la elección ni la puerta que contiene el premio. Eso es útil en el contexto del concurso, pero no es útil en absoluto a la hora de elaborar el código).

    Con este programa el experimento se repite 10 000 veces y arroja una probabilidad de 0.33, aproximadamente. Eso es 1/3. Pero si se modifica el código para cambiar de puerta una vez que se abre la que no tiene premio:

        # establecer el número de repeticiones del
        # experimento y configurar las puertas como antes:
    B <- 10000
    cambiar <- replicate(B, {
    	puertas <- as.character(1:3)
    	premios <- sample(c("coche","cabra","cabra"))
    	puerta_premiada <- puertas[premios == "coche"]
    	eleccion  <- sample(puertas, 1)
    	abrir <- sample(puertas[!puertas %in% c(eleccion, puerta_premiada)],1)
        # en lugar de plantarse, cambiar de puerta a una que
        # no sea ni la escogida (obvio, si no, no sería
        # cambiarse) ni la que ha sido abierta. Aquí sí que
        # tiene sentido la línea anterior.
    	cambiar <- puertas[!puertas %in% c(eleccion, abrir)]
        # comprobar si el resultado corresponde a la puerta premiada.
    	cambiar == puerta_premiada
    })
        # probabilidad de ganar plantándose
    mean(cambiar)
    

    Entonces la probabilidad después de repetir el experimento 10 000 veces es de 0.66, es decir, 2/3.

  3. Ganar partidos

    En el siguiente ejercicio, dos equipos se juegan un campeonato en una serie de 7 partidos, y ambos tienen un 50% de probabilidades de ganar cada encuentro. Si el equipo A ha perdido el primer partido de la serie, ¿cuál es la probabilidad de que ganen el campeonato? Ambos equipos tienen las mismas probabilidades de ganar, p = 0.5.

    Para conocer las probabilidades se puede hacer una simulación Monte Carlo, pero también se pueden calcular haciendo una lista de todas y cada una de las posibles soluciones. Para ello se puede hacer lo siguiente:

        # establecer el número de partidos que quedan por jugar
        # (los 7 del campeonato menos el primero, cuyo resultado
        # ya se conoce)
    n <- 7 - 1
        # establecer los resultados posibles de cada partido
        # individual
    resultados <- c(0, 1)
        # hacer una lista con los posibles resultados para 
        # el resto del campeonato repitiendo la lista hasta
        # conseguir un vector de longitud n
    l <- rep(list(resultados), length.out = n)
        # crear un conjunto de datos usando expand.grid
        # con todas las posibilidades combinatorias de los
        # elementos
    posibilidades <- expand.grid(l)
        # identificar cuáles son las filas (combinaciones de
        # partidos) que resultan en la victoria de A: como A perdió
        # el primer partido, los puntos necesarios siguen siendo 4.
    ganar <- rowSums(posibilidades) >= 4
        # comprobar la proporción de TRUE que hay en el vector 'ganar'
    mean(ganar)
    

    Por supuesto, también se puede hacer una simulación y repetirla un número grande de veces para comprobar el resultado. Para hacer una estimación Monte Carlo que dé el mismo resultado:

    B <- 10000 # establecer el número de repeticiones del experimento
    	   # definir el experimento a repetir
    A_ganan <- replicate(B, {
    	   # operaciones dentro del experimento:
    	   # hacer una serie aleatoria de 6 partidos:
    	   # de los resultados posibles, sacar 6 con reemplazo
        resultado <- sample(c(0, 1), 6, replace = TRUE)
        resultado <- sum(resultado) >= 4 #comprobar si hay victoria de A
        }
        )
    mean(A_ganan) # proporción de victorias en el vector
    

    Al final, la proporción de victorias anda en torno a 1/3, es decir, 0.33.

    Otro caso posible: ambos equipos no tienen las mismas posibilidades de ganar porque A es mejor que B y entonces p > 0.5. Entonces, aunque la simulación Monte Carlo es más o menos igual, hay que efectuar algunos cambios:

    1. No es posible hacer una repetición de los experimentos y ya, porque es necesario repetir los experimentos cambiando las probabilidades de victoria de ambos equipos.
    2. Dado el punto anterior, lo que hay que hacer es crear una función que haga una simulación Monte Carlo para cada par de probabilidades p, 1-p.
    3. Como las funciones personalizadas no se pueden aplicar a un vector elemento a elemento, hay que efectuar tal aplicación a través de la funcion sapply().

    Por lo tanto, el código es más o menos el siguiente:

    p  <- seq(0.5, 0.95, 0.025)
        # punto 2 de la lista anterior
    prob_B <- function(p){
        B <- 10000
        resultado <- replicate(B, {
    	B_gana <- sample(c(1,0), 7, replace = TRUE, prob = c(1-p, p))
    	sum(B_gana) >= 4
        })
        mean(resultado)
    }
        # punto 3 de la lista anterior
    Pr <- sapply(p, prob_B)
        # elaborar un gráfico para mostrar las probabilidades
        # de victoria para cada estimación de probabilidad
    png(file = "./incl/3-132f.png")
    plot(p, Pr)
    dev.off()
    

    3-132f.png

    Figure 17: Gráfico de probabilidades de ganar del equipo B

    Un ejemplo más para ver cómo funcionan las funciones personalizadas: en este caso se toma como probabilidad fija de que A gane el valor p = 0.75 y lo que se varía es la longitud de la serie de partidos. ¿Cuál es la probabilidad de ganar de B si su probabilidad de ganar un partido es 0.25 (es decir, A es 0.75) dadas series ficticias de 1, 3, 5 … 23, 25 partidos?

    prob_B <- function(N, p = 0.75){  # la función espera un valor N y toma
    				  # el valor de p como fijo
        B <- 10000
        resultado <- replicate(B, {
    	b_gana <- sample(c(1,0), N, replace = TRUE, prob = c(1-p, p))
    	sum(b_gana) >= (N+1)/2
    	})
          mean(resultado)
    }
    N <- seq(1, 25, 2)
    Pr <- sapply(N, prob_B)  # se aplica sobre todos los valores de N (el 
    			 # vector completo) la función prob_B()
    
    png(file = "plot.png")
    plot(N, Pr)
    dev.off()
    

    3-132g.png

    Figure 18: Gráfico del ejemplo anterior

4.2 Probabilidad continua

Cuando hay un conjunto de datos en el que cada elemento es una medida única, como en el conjunto heights, en el que las estaturas son todas diferentes, el asignar una probabilidad a cada medida como en el apartado anterior no es nada útil. Las probabilidades de que alguien esté en un grupo como los que vimos antes (ingesta de alcohol de 120 gramos diarios o más, por ejemplo) son útiles y significativas, pero cuando las medidas son únicas y exactas y no se repiten prácticamente en todo el conjunto de datos, como una altura de 1,72348 metros, las probabilidades de aparición son de 1/n 29 en realidad.

Igual que con las distribuciones, es mucho más útil definir una función que trabaje en intervalos en lugar de valores individuales. Repitiendo el código en aquel apartado (más o menos):

library(tidyverse)
library(dslabs)
data(heights)

    # los datos que se van a usar para la distribución
    # son las alturas de los varones.
datos <- heights %>%
	 filter(sex == "Male") %>%
	 pull(heights)

    # las medidas se toman desde el valor mínimo hasta el
    # valor máximo en 100 partes diferentes.
a <- seq(min(datos), max(datos), length = 100)

    # la función toma el valor de a y calcula la media de
    # los valores que son menores o iguales que a.
funcion_cdf <- function(a){
    mean(datos <= a)
}

    # la función sapply() aplica la funcion funcion_cdf()
    # a todos los valores de a (al vector completo)
valores_cdf <- sapply(a, funcion_cdf)

    # dibujar un gráfico simple
plot(a, valores_cdf)

Este código, que establece la distribución de valores del conjunto de datos, sirve para calcular las probabilides de aparición de una medida determinada. Por ejemplo, ¿cuál es la probabilidad de que la estatura de una persona tomada al azar sea superior a 70 pulgadas (las medidas de este conjunto de datos están en pulgadas)? Siguiendo el código de arriba, sería 1 - funcion_cdf(70)30.

En los conjuntos de datos que siguen una distribución normal, la probabilidad de que una medida esté por debajo de un cierto valor se calcula con la función pnorm() (ver este apartado). La densidad de la probabilidad se define como el área bajo la curva de densidad de la función que define el conjunto de datos (como la distribución normal). La densidad de probabilidad es otra fórmula matemática y es calculada por la función dnorm().31

4.2.1 Simulaciones Monte Carlo en distribuciones contínuas.

library(tidyverse)
library(dslabs)
data(heights)

x <- heights %>% filter(sex == "Male") %>%
    pull(height)

n <- length(x)
media <- mean(x)
sigma <- sd(x)

alturas_simuladas <- rnorm(n, media, sigma)

data.frame(alturas_simuladas = alturas_simuladas) %>%
    ggplot(aes(alturas_simuladas)) +
    geom_histogram(color = "black", binwidth = 2)

En este código lo que se hace es crear una simulación de las alturas reales usando la media y desviación típica del conjunto de datos reales. Del conjunto de datos heights se extrae el subconjunto de alturas de los varones. Luego se calculan la longitud del vector (cuántas medidas hay), la media y la desviación típica. Con estos datos se extraen tantas medidas de una distribución normal teórica como medidas hay en el conjunto de datos reales (n). La función rnorm() es lo que hace, extraer n medidas aleatorias de una distribución normal de media media y desviación típica sigma. Finalmente se crea un histograma para mostrar la forma de la distribución de datos.

Después viene la simulación propiamente dicha:

B <- 10000
alto <- replicate(B, {
    datos_simulados <- rnorm(800, media, sigma)
    max(datos_simulados)
})

mean(alto >= 7*12)

Se replica 10 000 veces una extracción del valor más alto de entre 800 generados aleatoriamente en una distribución normal de media media y desviación típica sigma. ¿Cuál es la probabilidad de que el más alto de un grupo de 800 sea más alto de 7 pies32? Dado que los valores que definen la distrbución normal han sido sacados de una muestra real, la simulación aproxima una muestra real, y el resultado, como era de esperar, es muy bajo: 0.019.

Volviendo a lo básico, las probabilidades en las distribuciones continuas se pueden calcular con las funciones dedicadas a las distribuciones continuas (valga la redundancia). La función norm() establece una distribución normal (teórica), por eso dentro de norm() encontramos (hay que fijarse en el prefijo de la función):

  • dnorm() calcula la densidad de probabilidad.
  • pnorm() calcula la probabilidad.
  • qnorm() calcula los cuantiles.
  • rnorm() calcula (extrae) valores aleatorios.

Hay otras distribuciones teóricas continuas, como student-t, chi-squared, beta, gamma o exponencial. R provee funciones para la mayoría de estas distribuciones teóricas.

4.3 Variables aleatorias, modelos de muestreo y el teorema del limite central

4.3.1 Variables aleatorias y modelos de muestreo

Una variable aleatoria es una variable numérica resultante de un proceso aleatorio. Un ejemple sencillo es el de la urna de bolas rojas y azules de ejemplos anteriores:

bolas <- rep(c("azul", "roja"), times = c(3, 2))
X <- ifelse(sample(bolas, 1) == "azul", 1, 0)

Es un código que ya se vio anteriormente. Equivale a escribir los colores de las canicas uno por uno. La función rep() replica un elemento (el vector en este caso) un cierto número de veces especificado en el parámetro times, en este caso tres veces el primer nombre y dos veces el segundo.

La segunda línea sí que es nueva. Es una sentencia condicional que asigna un 1 al elemento X si la extracción aleatoria de un elemento del vector bolas es azul y asigna un 0 en caso contrario (es decir, si es roja).

En ciencia de datos, a menudo hay que lidiar con datos afectados por el azar, como datos provenientes de una muestra aleatoria, datos afectados por errores de medida, o datos medidos en fuentes que son aleatorias en sí mismas. Es muy importante en este terreno ser capaz de medir la incertidumbre introducida por la aleatoriedad.

Muchos procedimientos de recolección de datos pueden ser modelados como extracciones de elementos de una urna, ya sean votantes, muestras aleatorias de la población de interés, etc. La urna contiene los valores posibles en esa población de interés de la que se extraerá la muestra. Los juegos de casino proporcionan buenos ejemplos de aplicaciones en el mundo real y son relativamente sencillos, por lo que es buena idea empezar por los juegos de casino.

  1. El ejemplo del casino y la ruleta

    Un casino me quiere contratar para saber si merece o no la pena instalar una ruleta. Si la probabilidad de pérdida es muy alta, no se instalará, pues debe dejar un margen de beneficio suficiente para pagar las facturas y las nóminas y que quede remanente para el enriquecimiento desmedido. Asumiremos aquí que van a jugar 1000 personas y que el único juego permitido es apostar a rojo o negro. Así que se trata de predecir cuánto dinero va a ganar o perder el casino.

    Vamos a definir una variable aleatoria, S, que representa los beneficios totales del casino. Para ello, primero hay que entender una ruleta (americana)33: esta consta de 36 números rojos o negros (18 de cada) y dos números verdes, y cada vez que se hace una tirada es como sacar los números de una urna con reemplazo. Así que una ruleta teórica se construye así:

    color <- rep(c("rojo", "negro", "verde"), times = c(18, 18, 2))
    

    Cada vez que sale rojo, la banca pierde 1€, y si no, gana 1€ (en realidad es más complicado, pero es para mantener el ejemplo sencillo). Una lista de 1000 jugadas se puede simular con la línea siguiente:

    color <- rep(c("rojo", "negro", "verde"), times = c(18, 18, 2))
    
    n <- 1000
    X <- sample(ifelse(color == "rojo", -1, 1), n, replace = TRUE)
    

    Esto es lo que se llama un modelo de muestreo, porque estamos modelando el comportamiento aleatorio de una ruleta muestreando la extracción de bolas (la ocurrencia de colores en las jugadas). Ahora que ya se han definido las 1000 jugadas y su resultado en una ruleta teórica, la variable aleatoria S que contiene las ganancias del casino es simplemente la suma del vector: S <- sum(X). Cada vez que se repite la ejecución del código el resultado varía porque S es una variable aleatoria, no tiende a ningún punto en particular.

  2. La distribución de probabilidad de la variable aleatoria

    Algo muy importante a tener en cuenta es la distribución de probabilidad de la variable aleatoria. Esta distribución de probabilidad nos ofrece la posibilidad de conocer la probabilidad de que el valor de la variable se encuentre en un determinado intervalo. Por ejemplo, para conocer la probabilidad de perder dinero (el casino), se puede medir la probabilidad de S < 0.

    Si se puede definir una función de distribución acumulativa o CDF (ver aquí) entonces se puede responder cualquier pregunta sobre las probabilidades contenidas en S.

    F(a) = Pr(S <= a)

    Haciendo una simulación Monte Carlo de un número grande de secuencias de 1000 jugadas individuales se podrá hacer una estimación de cómo se distribuyen estos resultados.

    n <- 1000
    B <- 10000
    S <- replicate(B, {
        X <- sample(c(-1, 1), n, replace = TRUE,
    		prob = c(9/19, 10/19))
        sum(X)
    }
    )
    

    Un sencillo histograma hist(S) muestra la distribución de las ganancias:

    3-312b.png

    Y calculando la media de las veces que el casino pierde dinero, se conoce la probabilidad exacta: mean(S < 0), exactamente 0.0458 o 4.58%.

    Se puede ver que la distribución es prácticamente normal (un diagrama Q-Q lo puede confirmar). Por lo tanto, si se ajusta a lo que es una distribución normal, so pueden conocer los detalles con solo la media y la desviación típica. Ahora se puede crear un gráfico para ver la distribución y la densidad de probabilidad.

    library(tidyverse)
    s <- seq(min(S), max(S), length = 100)    # dividir el rango de valores de S en 100 partes
    normal_density <- data.frame(s = s, f = dnorm(s, mean(S), sd(S))) 
    data.frame (S = S) %>%    # make data frame of S for histogram
        ggplot(aes(S, ..density..)) +
        geom_histogram(color = "black", binwidth = 10) +
        ylab("Probabilidad") +
        geom_line(data = normal_density, mapping = aes(s, f), color = "blue")
    

    Se divide el rango de valores de S (del mínimo al máximo) en 100 partes. Luego se crea un conjunto de datos con estos valores de S y con la distribución de densidad normal para S, es decir, basada en sus valores de media y desviación típica. Luego, con el conjunto de datos se crea el gráfico.

    NOTA IMPORTANTE SOBRE DISTRIBUCIONES Y PROBABILIDAD Si tenemos una lista de números x _{1}, x _{2}, … , x _{n} la distribución es una función que se define como la proporción de números de la lista que son menores que una cantidad determinada F(a) = x <= a. Sobre esta lista se pueden hacer operaciones para calcular valores muy útiles cuando la distribución se aproxima a una distribución normal, como la media o la desviación típica.

    Una variable aleatoria X tiene una función de distribución, pero no se necesita una lista de números porque es un concepto teórico. Se define la función F(a) que responde a la pregunta de cuál es la probabilidad de que X sea igual o menor que a?.

    Sin embargo, si X se define mediante el proceso de sacar números de una urna, entonces que hay una lista: la lista de números dentro de la urna. Entonces, su distribución (la distribución de la variable aleatoria X) es la de la lista de números de la urna. La media de la lista de números es el valor esperado de la variable aleatoria X y la desviación típica de la lista de números es el error estándar de X.

    También puede hacerse una simulación Monte Carlo de una variable aleatoria. Aquí también hay una lista de números que aproximará más o menos exactamente los valores esperados de X. Cuanto más grande el número de simulaciones, mejor será la aproximación. Nuevamente, el valor esperado de X es la media de la lista de números y el error típico es la desviación típica de la lista.

    Es decir, según yo lo entiendo, es que una variable aleatoria tiene un valor esperado (que es la media de los valores x de la lista de números que aparecen en la urna o en la simulación) y un error típico. O sea, que la variable vale Y ± y siendo Y el valor esperado e y el error típico, y la distribución de probabilidad de que así sea es la distribución de la lista de números.

    De manera que se pueden definir una serie de conceptos matemáticos muy útiles a partir de estos principios, como que el valor esperado de una variable aleatoria que se defina por una extracción independiente de un valor34 es la media de los valores existentes. Por ejemplo, en la ruleta que se modeló más arriba para simular una ruleta y apostar al rojo, el valor esperado sería

    E(X) = (20 + (-18)) / 38

    Es decir, el valor esperado (E) de la variable aleatoria X es 20 euros positivos resultado de los números rojos y verdes más 18 euros negativos resultado de sacar los números negros, dividido entre 38 números en total (la media, vamos), lo cual arroja un resultado de 5.2 céntimos (exactamente 0.05263158). Parece en principio un poco extraño que el valor esperado sea 0.05 cuando los valores existentes son 1 y -1. Si hubiera 19 euros positivos y 19 euros negativos, el valor esperado sería 0, pero como los valores no están equilibrados, la media es diferente. En realidad significa que de cada euro jugado, el casino gana, de media, 5 céntimos, lo cual sería la media de las extracciones sucesivas: cuantas más extracciones, más cercano el valor a la media teórica.

    La fórmula general para la media (el valor esperado) de una variable aleatoria con dos posibles resultados, a y b cada uno de los cuales tiene una probabilidad de aparición p y (1 - p) respectivamente es

    E(X) = ap + b (1 - p)

4.3.2 El teorema del límite central

Cuando la muestra es suficientemente grande, la suma de las comprobaciones aleatorias de la muestra es aproximadamente normal.

En fin, sabemos que si una lista de números es aproximadamente normal entonces puede ser descrita con solo dos números: la media y la desviación típica. También sabemos que lo mismo se aplica a la distribución de probabilidad de una variable aleatoria. En esta última, los números que la describen son el valor esperado y el error típico.

El valor esperado de una variable aleatoria X es, tal y como se vio más arriba:

E(X) = μ

Esta definición del valor esperado es útil para aproximar la distribución de probabilidad de las sumas de los valores extraídos en cada extracción o juego, lo que a su vez es útil para describir las distribuciones de medias y proporciones.

Lo primero es considerar que el valor esperado de la suma de las extracciones es el número de extracciones efectuadas multiplicado por la media de valores en la muestra (la urna o la ruleta, por ejemplo). Así que en el caso de la ruleta y el casino, dado que el valor esperado de la variable aleatoria es 5 céntimos, si 1,000 personas juegan a la ruleta, el casino ganará en promedio 50€.

Al mismo tiempo, se puede esperar (se debe esperar) que cada uno de los valores individuales, cada una de las observaciones hechas difiera de esta media. De hecho, en el ejemplo de la ruleta no cabe otra cosa, pues aunque la media es 0.05, los valores son o 1 o -1, es imposible que aparezca un 0.05, no existe tal posibilidad. ¿Cuál es la probabilidad de aparación de los distintos valores?

El error típico (SE) nos da una idea de cómo son las variaciones alrededor del valor esperado (que en el caso de la ruleta, repito, no existe como observación individual). Si las observaciones son independientes (importante), entonces el error típico es la raíz cuadrada del número de obserbaciones por la desviación típica de los valores de la urna:

SE[X] = √observaciones x σ

Y para conocer la desviación típica, hay que conocer las probabilidades de aparición de cada uno de los valores. Es decir, si en la muestra hay dos valores, a y b con probabilidades de aparición p y (1-p) respectivamente, entonces la desviación típica será el valor absoluto de a - b multiplicado por la raíz cuadrada de p (1-p), o lo que es lo mismo

a - b * √p(1-p)

Lo que arroja un resultado en el ejemplo de la ruleta de |1-(-1)|*√(10/19 - 9/19) = 2√90/19 = 0.99 lo que es prácticamente 1: o sale 1 o sale -1. Por lo que cada observación será el valor esperado de la variable (0.05) ± 1 con un resultado ligeramente favorecedor del 1 frente al -1.

Siguiendo la fórmula anterior que define el valor de SE, se puede calcular que √1000 · 2 · √(90)/19 es 31.58€. Resultado final, en 1,000 jugadas el casino espera una media de 50€ ±31.6€.

Ahora, sabiendo que el teorema del límite central dice que la suma de los valores es aproximadamente normal nos podemos saltar las simulaciones Monte Carlo y usar la aproximación normal para conocer la probabilidad de que el casino pierda dinero.

n <- 1000
mu <- n * (20-18)/38
se <- sqrt(n) * 2 * sqrt(90)/19

pnorm(0, mu , se)

La probabilidad de que el valor de la función normal de media μ y error típico σ sea igual o menor que 0 es de menos del 5%.

4.4 Propiedades estadísticas del promedio o media

4.4.1 Propiedad nº 1

El valor esperado de la suma de varias variables aleatorias es la suma de los valores esperados de las mismas variables.

Este es el caso general, pero en un aspecto más particular se puede inferir que si las variables aleatorias son observaciones independientes de la misma variable, entonces se cumple la norma explicada más arriba, que la suma de las observaciones es el número de observaciones por el valor esperado de cada observación (siempre el mismo, pues son observaciones independientes de la misma urna). Es decir, n * E[X].

4.4.2 Propiedad nº 2

Una variable multiplicada por una constante no aleatoria es igual a la constante por el valor esperado original de la variable.

E[aX] = a * E[X]

La manera más intuitiva de entenderlo es el cambio de unidades. Si el valor esperado es 0.05€, también puede multiplicarse por 100 para medirlo en céntimos de euro sin que haya ningún cambio.

Una consecuencia de las dos propiedades anteriores es que el valor esperado de la media de observaciones independientes de la misma muestra (la urna) es el valor esperado de la muestra, o sea μ.

E[(X _{1} + X _{2} + · · · + X _{n} )/n] = E[X _{1} + X _{2} + · · · + X _{n}]/n = nµ/n = µ

4.4.3 Propiedad nº 3

El cuadrado del error típico de la suma de variables aleatorias independientes es la suma de los cuadrados de los errores típicos de las variables una por una.

El cuadrado del error típico es conocido habitualmente como varianza o variancia, σ2.

Es lo mismo que decir que el error típico de la suma de variables aleatorias independientes es la raíz cuadrada de la suma de los cuadrados de los errores típicos de las variables.

4.4.4 Propiedad nº 4

El error típico de una variable aleatoria multiplicada por una constante es igual a la constante multiplicada por el error típico de la variable.

SE[aX] = a · SE[X]

Una consecuencia de las propiedades 3 y 4 es que el error típico de la media de extracciones independientes de una urna es la desviación típica de los elementos de la urna dividida por la raíz cuadrada de n (el número de extracciones).

4.4.5 Propiedad nº 5

Si X es una variable aleatoria normalmente distribuida y a y b son dos constantes, entonces aX + b también es normalmente distribuida.

Todo lo que se consigue es cambiar las unidades (como en propiedades anteriores) al multiplicar, y desplazar el centro de la distribución normal al sumar un valor. Pero la forma de la distribución sigue siendo la misma.

4.5 La ley de los grandes números

Una implicación muy importante de las propiedades descritas más arriba es que a medida que crece el número de observaciones, el error típico desciendo más y más, hasta que finalmente la media de las observaciones converge con la media de la muestra.

Esta es una norma que se suele malinterpretar al considerar que si salen 5 caras en 5 tiradas de una moneda, es más probable que la sexta sea cruz para equilibrar la media esperada del 50%. Pero esta ley solo se aplica a números grandes, no pequeños. Después de 5 caras, la probabilidad de cara en la siguiente tirada sigue siendo del 50%.

El Teorema del Límite Central funciona con un número grande de observaciones, pero «grande» es un término relativo. En ocasiones en las que la probabilidad del resultado es especialmente alta, 40 observaciones pueden ser suficientes, incluso 10. En otros casos, cientos de miles de observaciones pueden no ser suficientes. Es el caso de la lotería, a pesar de haber cientos de miles de boletos, solo 1 o 2 salen premiados. Para aproximar normalmente esta distribución serían necesarios muchos millones de observaciones, y aún así… Para casos como este, sería más apropiada una distribución de Poisson.

4.6 Tasas de interés

Los cálculos anteriores también pueden ser usados por los bancos para calcular la tasa de interés que se aplica a un préstamo.

Si el banco sabe que aproximadamente un 2% de los deudores falla en pagar su deuda, el préstamo sin intereses hará que incurra en pérdidas, pues el banco no sabe exactamente quiénes son los que pagarán y los que no. La aportación extra que proporcionan los intereses hará que no pierda dinero, pueda pagar las nóminas y tener un beneficio. Pero si el interés es demasiado alto los clientes pueden irse a otro banco.

Por ejemplo, el banco presta 180 000€ a 1 000 personas este año. Supongamos, por simplificar, que se pierden 200 000€ por ejecuciones hipotecarias y que esta cifra incluye todos los gastos de operaciones. Dados estos supuestos, el escenario se puede modelar tal que así:

n <- 1000
perdidas_por_ejecucion <- -200000
p <- 0.02

impagos <- sample(c(0, 1), n, prob = c(1-p, p), replace = TRUE)

sum(impagos * perdidas_por_ejecucion)

Como cada vez que no se paga un préstamo el banco pierde 200 000€, multiplicar estas pérdidas por la lista de ceros y unos que resulta de modelar una lista de morosos (los unos) entre los que pagan (los ceros) y sumar el resultado nos da una idea de las pérdidas globales.

Cada vez que se ejecuta el código, como está basado en probabilidades y no números ciertos, el resultado será distinto, pero siempre habrá pérdidas (gracias a ese 2% del que se habló antes). Una simulación Monte Carlo nos ofrece una visión de la distribución de la variable:

library(tidyverse)

n <- 1000
perdidas_por_ejecucion <- -200000
p <- 0.02

B <- 10000

perdidas <- replicate(B, {
   impagos <- sample(c(0, 1), n, prob = c(1-p, p), replace = TRUE)
   sum(impagos * perdidas_por_ejecucion)
})

data.frame(perdidas_en_millones = perdidas / 10^6) %>%
   ggplot(aes(perdidas_en_millones)) + 
       geom_histogram(binwidth = 0.6, col = "black")

En realidad, esto no es necesario, pues el teorema del límite central dice que la suma de una variable aleatoria está normalmente distrubuida y puede ser descrita con los valores de media y desviación típica siguiendo las fórmulas que ya conocemos:

mu <- n * (p * perdidas_por_ejecucion + (1-p) * 0)
sigma <- sqrt(n) * abs(perdidas_por_ejecucion - 0) * sqrt(p * (1-p))

En fin, que se pierde dinero sí o sí.

Lo que el banco desea es hallar un valor x aplicable como beneficio opuesto a las posibilidades de pérdida (es decir, a los préstamos que sí se devuelven) que satisfaga la siguiente ecuación:

l · p + x (1 - p) = 0, donde l es el valor de las pérdidas atribuibles a cada ejecución.

Es decir, el valor mínimo (por eso lo igualamos a 0, para asegurar que no hay pérdidas) con el que no hay pérdidas, pero tampoco ganancias. Sustituyendo en la ecuación anterior, el valor de x es

x = -(l · p) / (1 - p)

o lo que es lo mismo

x <- -perdidas_por_ejecucion * p / (1-p)

El resultado, 4081.633, es aproximadamente un 2% de estas pérdidas por ejecución (exactamente un 2.04%). Aplicando un interés del 2%, la suma de pérdidas y ganancias será aproximadamente cero.

Algo más complejo es asegurarnos de que la probabilidad de perder dinero es del 1% solamente: Pr(S < 0) = 0.01. Sabemos que S es aproximadamente normal y que el valor esperado de S es

E[S] = {l · p + x (1 - p)} · n

y su error típico es

SE[S] = | x - l | · √p (1 - p) · √n

Ahora se usa un truco matemático que consiste en sumar y dividir ambos términos en la primera ecuación para que no cambie la igualdad. Se trata de convertir

Pr(S < 0) = 0.01

en

Pr((S - E[S] / SE[S]) < (-E[S]/SE[S]))

Esto tiene un par de ventajas, y es que el término de la izquierda es una variable aleatoria normal que llamaremos Z con valor esperado 0 y error típico 135. La otra ventaja es que si lo anterior es cierto, entonces el término de la derecha es igual a qnorm(0.01). Esto es así porque una variable normal cuyo valor esperado es igual o menor que 0.01 es la definición de distribución normal. Después de revisar el uso de qnorm() y pnorm() más arriba, veo que esto es porque qnorm(0.01) es el valor de Z para el que pnorm() es 0.01.

De manera que sustituyendo E[S] y SE[S] en la ecuación anterior con los valores un poco más arriba, y despejando la x se obtiene que36

l <- perdidas_por_ejecucion
z <- qnorm(0.01)
x <- -l*( n*p - z*sqrt(n*p*(1-p)))/ ( n*(1-p) + z*sqrt(n*p*(1-p)))
x/180000    # tasa de interes
perdidas_por_ejecucion*p + x*(1-p)    # valor esperado del beneficio por
				      #  prestamo
n*(perdidas_por_ejecucion*p + x*(1-p)) # valor esperado del beneficio global
				       # en n préstamos

El resultado es 6249, o sea un 3.5%. El valor esperado anual (en los 1 000 préstamos) es de n · ( l · p + x · (1 - p)), es decir 2 124 198€.

Una simulación Monte Carlo confirma esto después de definir las variables con el código anterior:

B <- 100000
beneficio <- replicate(B, {
    prestamos <- sample( c(x, perdidas_por_ejecucion), n, 
			prob=c(1-p, p), replace = TRUE) 
    sum(prestamos)
})
mean(beneficio)    # expected value of the profit over n loans
mean(beneficio < 0)    # probability of losing money

4.7 Tasas de interés II

En realidad, una probabilidad «estándar» para todos y cada uno de los préstamos o pólizas de seguro es poco realista, lo lógico es que haya variaciones para cada uno de los contratadores o prestatarios. De ahí que el código se varíe un poco para incluir una ligera variación en la tasa de probabilidad que se aplica, una variación para cada uno de los préstamos.

  p <- 0.015
  l <- -150000  # pérdidas por cada impago
  n <- 1000
  B <- 10000
  z <- qnorm(0.05) # posibilidad de pérdidas del 5%
  x <- -l*( n*p - z*sqrt(n*p*(1-p))) / ( n*(1-p) + z*sqrt(n*p*(1-p))
     # x es el valor de la póliza necesario para cumplir con los
     # requisitos anteriores

     # ¿Cuáles son los beneficios si las probabilidades de defunción (p)
     # varían individualmente ±0.01?
  beneficio <- replicate(B, {
      new_p <- p + sample(seq(-0.01, 0.01, length = 100), 1)
      observaciones <- sample(c(l, x), n,
			      prob = c(new_p, 1 - new_p),
			      replace = TRUE)
       sum(observaciones)
})

     # Resultados
mean(beneficio)  # beneficio esperado
mean(beneficio < 0) # probabilidad de perder dinero
mean(beneficio < -1e6) # prob. de perder más de 1 millón

5 Inferencia estadística

6 Otras funciones interesantes

6.1 letters

Esta función devuelve las letras del alfabeto. En realidad no es una función (creo), sino una variable predefinida. Así, letters[2], es decir, la segunda posición del vector letters devuelve "b", y letters[3:5] devuelve "b" "c" "d" "e".

6.2 duplicated

Esta función crea un vector lógico de la misma longitud que el original en el que aparecen TRUE en la posición en la que un elemento del vector original aparece una segunda vez (o posteriores). Por lo tanto, si el vector original es 1 2 3 4 3 1, el vector resultado de aplicar la función será FALSE, FALSE, FALSE, FALSE, TRUE, TRUE.

6.3 any

Es TRUE si en el vector lógico pasado como parámetro hay al menos un valor TRUE. Será FALSE si no hay ninguno.

6.4 sapply

Dado que algunas funciones no se aplican a un vector elemento a elemento, como por ejemplo las creadas ad hoc, hace falta una función que fuerce esta aplicación. La función sapply(x, f) hace que la función f se aplique al vector x de esta manera elemento por elemento.

Footnotes:

1

En realidad esto no es necesario, cualquier extensión vale, pero los sistemas operativos entenderán la extensión .R para colocar el icono correspondiente, abrir el IDE con un doble clic, etc.

2

El orden en el que aparecen los factores (las cuatro categorías o, en este caso, regiones) no respeta el orden de aparición dentro del vector, sino que aparecen en orden alfabético. Se pueden reordenar con la función reorder(), pero se verá en la segunda parte: visualización de datos.

3

Porque el vector tiene asociados nombres a sus elementos (con la función names).

4

Creo que coerción está bien, pero cuando se trata de un adjetivo (como en El vector coercionado es el último de la lista) se entiende mejor forzado.

5

Puede que haya que instalarlo: install.packages("dplyr") y cargarlo luego con library(dplyr).

6

El recíproco de x es 1/x o x-1.

7

Por supuesto, i es una convención, se puede usar cualquier variable.

8

La funcíon table muestra la frecuencia de aparición de cada valor individual de un conjunto de datos.

9

¿Por qué? Buena pregunta.

10

Hay que estudiar esto más a fondo. No sé lo que es.

11

También valdría print(p).

12

Todas las funciones siguen el mismo esquema: geom, un guion bajo y un nombre alusivo al tipo de gráfico que se quiere conseguir.

13

Comprobación de los argumentos posibles de una función: args(ggplot).

14

En una consola de texto, cd . no cambia de directorio, pues el punto representa el directorio actual. Por eso, una orden como cp ./carpeta/archivo.ext otro/directorio/diferente/ busca (para copiarlo) el archivo archivo.ext dentro de la carpeta carpeta, que está dentro del directorio actual.

15

La función topn() ha sido superada por slicemin() y slicemax() (y otras). Aunque sigue funcionando, no recibe más actualizaciones que las críticas y desde luego ninguna funcionalidad nueva. Todo el desarrollo de código se ha trasladado a slice.

16

En cualquier caso, lo de 'tallado' es idea mía. Es posible que la costumbre en la ciencia de datos en español sea darle otro nombre o que se mantenga el original inglés faceting. En cuanto me entere, lo corrijo.

17

Es útil consultar la ayuda de esta función para conocer las opciones que ofrece: help(themes).

18

Creo que las mejoras se deben a la fuerza en los datos de las mejoras en China, Korea y, probablemente, Japón (y otros países del Extremo Oriente), el resto están jodidos e incluso peor que antes. Mientras tanto, el grupo Occidente ha pasado de 32$ diarios per cápita a 64. No está mal.

19

El paquete RColorBrewer proporciona paletas de colores adecuadas para cada una de estas ocasiones. Más información sobreestas paletas en el libro de texto, en la sección Encoding a third variable. El siguiente código muestra las paletas de colores:

library(RColorBrewer)
display.brewer.all(type = "seq")

Usar type = "div" muestra los colores para patrones divergentes.

20

Recuérdese que dos dígitos significativos son un solo decimal.

21

Las probabilidades de que Obama ganase las elecciones eran 80% en 2008 y 90% en 2012, y ganó. Para Hillary Clinton la probabilidad era del 71% y perdió.

22

En este caso tampoco pasaría nada, son cinco bolas nada más, pero si fueran 500…

23

Aquí se trata la baraja francesa: As, 2-10, Sota, Reina, Rey (13 cartas) y cuatro palos: Picas, Diamantes, Tréboles y Corazones. Total, 52 cartas. Se pueden adaptar los ejemplos a la baraja española, pero no lo voy a hacer.

24

Si no se especifica un «título» para las variables (es decir, expand.grid(pant, cams), la tabla resultante mostrará Var1 Var2 como cabeceras de las columnas (Variable 1 y Variable 2).

25

En realidad, no son todos los posibles, pues no se repiten números y en los números de teléfono reales sí que está permitida la repetición.

26

En Blackjack, para sacar 21, que es el objetivo, hay que juntar un as y una figura (eso da 21 puntos directamente y es victoria segura). Aquí no importa sacar as y figura que figura y as. Tanto una combinación como la otra son válidas y no se requiere contarlas dos veces.

27

Una solución sería crear un bucle for para aplicar la función sobre n, pero siempre es preferible aplicar operaciones sobre vectores completos.

28

Un tal Monty Hall fue el más famoso de los presentadores de este concurso que se llamaba Trato Hecho Let's Make a Deal.

29

Donde n es el número de elementos del conjunto de datos.

30

La función sobre el valor 70 calcula la probabilidad de que las medidas sean menores de 70 y por lo tanto hay que restar de 1.

31

Ver el apartado 13.11.2 del libro de texto.

32

7 pies son 7*12 pulgadas, o sea, 84 pulgadas. En cristiano, 2,13m. Los valores del conjunto de datos están en pulgadas, de ahí esta aclaración.

33

Una ruleta americana tiene 18 números rojos, 18 negros y 2 verdes (el 0 y el 00). Da mayor márgen de beneficio a la banca. La ruleta europea solo tiene un número verde, el 0.

34

Como la urna o la ruleta, por ejemplo, donde hay unos valores fijos «guardados» y se saca o se decide uno por vez.

35

Esto se me escapa. Me lo creo, pero tengo que investigarlo.

36

Es mejor seguir este razonamiento algebraico en el libro, página 284, sección 14.11, porque es una locura escribirlo aquí. Y mi álgebra está demasiado oxidada.

Author: JAC

Created: 2021-09-14 mar 21:23

Validate