Violin plot en R

Los gráficos de violín, o violin plots, son una alternativa a los diagramas de caja que resuelven los problemas relacionados con la visualización de la distribución subyacente de los datos, ya que estos gráficos muestran una estimación de la densidad tipo núcleo de las observaciones. En este tutorial mostraremos cómo crear un un violin plot en R base a partir de un vector y de data frames, cómo agregar puntos que representan la media y dividir (partir) los gráficos de violín.

Gráfico de violín de un vector

Para crear un violin plot en R a partir de un vector, debes pasar el vector a la función vioplot del paquete del mismo nombre. Considera, por ejemplo, el siguiente vector:

x <- c(6, 9, 0, 19, -1, 8, 12, 5, 3, 7,
       2, 4, 3, -8, -9, 8, 4, 12, 5, 14)

Usando el vector de ejemplo, puedes crear un simple diagrama de violín en R escribiendo:

# install.packages("vioplot")
library("vioplot")

vioplot(x)

Por defecto, la función vioplot crea un violín vertical, pero si estableces el argumento horizontal como TRUE, puedes crear un violin plot horizontal.

vioplot(x, horizontal = TRUE)

Si quieres personalizar el violin plot, existen varios argumentos para controlar los parámetros gráficos:

vioplot(x,
        col = 2,               # Color del área
        rectCol = "red",       # Color del rectángulo
        lineCol = "white",     # Color de la línea
        colMed = "green",      # Color del símbolo pch
        border = "black",      # Color del borde del violín
        pchMed = 16,           # Símbolo pch para la mediana
        plotCentre = "points") # Si "line", dibuja una línea en la mediana

Además, puedes añadir las observaciones al gráfico con la función stripchart como sigue:

stripchart(x, method = "jitter", col = "blue",
           vertical = TRUE, pch = 19, add = TRUE)

Ten en cuenta que si tienes un violin plot horizontal tendrás que especificar vertical = FALSE en la función anterior.

Cabe destacar que también es posible dibujar un violin plot sin tener en cuenta los datos atípicos. Para ese propósito, puedes asignar a una variable la salida de la función boxplot y luego devolver los valores del vector original que no son valores atípicos.

box <- boxplot(x)

# Nos quedamos con los datos no atípicos
x <- x[!(x %in% box$out)]

vioplot(x)

También puedes establecer el argumento ylog como TRUE en caso de que quieras que el eje Y esté en escala logarítmica. Ten en cuenta que solo se aceptan datos positivos.

par(mfrow = c(1, 2))

vioplot(1:10)
vioplot(1:10, ylog = TRUE)

par(mfrow = c(1, 1))

Histograma y violin plot

Por último, ten en cuenta que puedes dibujar un violin plot sobre un histograma. Considera, por ejemplo, que la distribución subyacente de tus datos presenta multimodalidad. En este caso, un diagrama de caja y bigotes no representa esta condición, mientras que un gráfico de violín sí lo hará. La siguiente representación gráfica te ayudará a entender mejor por qué un violin plot es útil:

set.seed(1)

# Datos multimodales
n <- 10000
ii <- rbinom(n, 1, 0.5)
datos <- rnorm(n, mean = 130, sd = 10) * ii +
         rnorm(n, mean = 80, sd = 5) * (1 - ii)

# Histograma
hist(datos, probability = TRUE, col = "grey", axes = FALSE,
     main = "", xlab = "",  ylab = "")

# Eje X
axis(1)

# Densidad
lines(density(datos), lwd = 2, col = "red")

# Añadimos el violin plot
par(new = TRUE)
vioplot(datos, horizontal = TRUE, yaxt = "n", axes = FALSE,
        col = rgb(0, 1, 1, alpha = 0.15))

Violin plot por grupo en R

Por una parte, si tienes un data frame donde una variable contiene grupos puedes crear un violin plot con una fórmula, especificando la variable numérica contra el factor. A continuación mostramos un ejemplo usando el conjunto de datos chickwts de R base, que contiene las variables weight (peso) y feed (tipo de alimentación).

tail(chickwts) # Últimas filas

    weight  feed
66    352  casein
67    359  casein
68    216  casein
69    222  casein
70    283  casein
71    332  casein

Ahora puedes especificar la fórmula en el primer argumento, los colores y demás parámetros gráficos que quieras:

data <- chickwts

vioplot(data$weight ~ data$feed, col = 2:length(levels(data$feed)),
        xlab = "Feed", ylab = "Weight")

También puedes añadir las observaciones al violin plot anterior con la función stripchart escribiendo lo que sigue:

stripchart(data$weight ~ data$feed, vertical = TRUE, method = "jitter",
           pch = 19, add = TRUE, col = 3:8)

Por otra parte, si tu conjunto de datos contiene columnas numéricas que representan algunas variables puedes crear el violin plot en R directamente del data frame. Como ejemplo, usaremos el conjunto de datos trees.

tail(trees) # Últimas filas

   Girth  Height  Volume
26  17.3    81     55.4
27  17.5    82     55.7
28  17.9    80     58.3
29  18.0    80     51.5
30  18.0    80     51.0
31  20.6    87     77.0

Si pasas el data frame a la función vioplot puedes crear el gráfico. Ten en cuenta que si apilas este data frame con la función stack, podrás especificar una fórmula igual que en el ejemplo anterior.

data <- trees

vioplot(data, col = 2:4, border = 2:4)

# Equivalente a:
datos_apilados <- stack(trees)
vioplot(datos_apilados$values ~ datos_apilados$ind, col = 2:4,
        border = 2:4)

Reordenar un violin plot

Los gráficos de violin en R se ordenan por defecto por el orden de los niveles de la variable categórica. Revisa el gráfico que creamos antes con el conjunto de datos chickwts y comprobarás que las variables siguen el siguiente orden:

levels(chickwts$feed)

 "casein"  "horsebean"  "linseed"  "meatmeal"  "soybean"  "sunflower"

Sin embargo, puedes sobrescribir este comportamiento reordenando la variable categórica por cualquier característica de los datos con la función reoder. En el siguiente ejemplo vamos a utilizar la mediana, pero puedes elegir la función que prefieras.

par(mfrow = c(1, 2))

data <- chickwts

#----------------
# De menos a más
#----------------

medianas <- reorder(data$feed, data$weight, median)
# medianas <- with(data, reorder(feed, weight, median)) # Equivalente

vioplot(data$weight ~ medianas, col = 2:(length(levels(data$feed)) + 1),
        xlab = "", ylab = "Weight", las = 2)

#----------------
# De más a maneos
#----------------

medianas <- reorder(data$feed, -data$weight, median)
# medianas <- with(data, reorder(feed, -weight, median)) # Equivalente

vioplot(data$weight ~ medianas, col = 2:(length(levels(data$feed)) + 1),
        xlab = "", ylab = "Weight", las = 2)

par(mfrow = c(1, 1))

Añadir media a un diagrama de violín

La función vioplot muestra la mediana de los datos, pero si la distribución no es simétrica la media y la mediana pueden ser muy diferentes. Por tanto, puede resultar de interés añadir la media, u otra métrica, a un violin plot en R base con la función points. Ten en cuenta que los pasos a seguir son diferentes si estás trazando un violin plot horizontal o vertical, o si estás creando gráficos individuales o por grupo.

Por una parte, para mostrar la media en un gráfico de violín creado con un solo vector, puedes escribir:

par(mfrow = c(1, 2))

# Datos exponenciales
set.seed(5)
x <- rexp(20)

#-------------------
# vioplot vertical
#-------------------

vioplot(x, col = 4)

# Añadimos la media
points(mean(x), pch = 19, col = "green", cex = 1.1)

#-------------------
# vioplot horizontal
#-------------------

vioplot(x, col = 4, horizontal = TRUE)

# Añadimos la media
points(mean(x), 1, pch = 19, col = "green", cex = 1.1)

legend("topright", pch = c(21, 19), col = c("black", "green"),
       bg = "white", legend = c("Mediana", "Media"), cex = 1.1)

par(mfrow = c(1, 1))

Por otra parte, si tienes un violin plot por grupos, para agregar la media puedes ejecutar lo que sigue:

par(mfrow = c(1, 2))

set.seed(5)
df <- data.frame(x = rexp(20), y = rexp(20), z = rexp(20))

#---------------------------
# vioplot vertical por grupo
#---------------------------

vioplot(df, col = 2:4)

# Añadimos las medias
means <- apply(df, 2, mean)
means <- colMeans(df) # Equivalente (más eficiente)

points(means, pch = 19, col = "green", cex = 1.1)

legend("top", pch = c(21, 19), col = c("black", "green"),
       bg = "white", legend = c("Mediana", "Media"), cex = 1.1)

#-----------------------------
# vioplot horizontal por grupo
#-----------------------------

vioplot(df, col = 2:4,
        horizontal = TRUE)

# Añadimos las medias
means <- apply(df, 2, mean)
means <- colMeans(df) # Equivalente (más eficiente)

points(means, 1:ncol(df), pch = 19, col = "green", cex = 1.1)

par(mfrow = c(1, 1))

Dividir violin plots en R

Para finalizar, cabe mencionar que puedes dividir un violin plot en R. Considera, por ejemplo, que has dividido el conjunto de datos trees en dos grupos, representando árboles altos y bajos dependiendo de su altura. En este caso puedes hacer uso de los argumentos side y add de la siguiente manera:

data <- trees

altos <- trees[trees$Height >= 76, ]
bajos <- trees[trees$Height < 76, ]

vioplot(altos, side = "left", plotCentre = "line", col = 2)
vioplot(bajos, side = "right", plotCentre = "line", col = 3, add = TRUE)

legend("topleft", legend = c("Altos", "Bajos"), fill = c(2, 3), cex = 1.25)