Gráfico de densidad en R

La función de densidad de probabilidad de un vector x, que se suele denotar como f(x) describe la probabilidad de que la variable tome un determinado valor. La función de densidad de probabilidad empírica es una versión suavizada del histograma. Esta suavización también se conoce como estimador de Parzen-Rosenblatt o estimador kernel. Puedes hacer un gráfico de densidad en R de manera muy sencilla tal y como mostraremos en este tutorial. Al final de la lectura sabrás cómo dibujar diagramas de densidad tanto en R como en RStudio.

Función de densidad empírica en R

Para crear un gráfico de densidad en R, puedes pasar a la función plot un objeto creado con la función density de R, que dibujará una curva de densidad en una nueva ventana. Si lo prefieres, también puedes superponer la línea de densidad sobre un histograma con la función lines.

set.seed(1234)

# Generamos datos
x <- rnorm(500)

par(mfrow = c(1, 2))

# Creamos un histograma
hist(x, freq = FALSE, main = "Histograma y densidad",
     ylab = "Densidad")

# Calculamos la densidad
dx <- density(x)

# Añadimos la línea de densidad
lines(dx, lwd = 2, col = "red")

# Curva de densidad sin histograma
plot(dx, lwd = 2, col = "red",
     main = "Densidad")

# Añadimos los datos con riudo en el eje X
rug(jitter(x))

El resultado es la densidad de probabilidad empírica suavizada. Una alternativa para crear una función de densidad en R es la función epdfPlot del paquete EnvStats. Con esta función puedes pasar el vector numérico directamente como parámetro.

# Alternativa equivalente con el paquete EnvStats
# install.packages("EnvStats")

library(EnvStats)

epdfPlot(x, epdf.col = "red")

Selección de ventana para densidades tipo núcleo

Cuando dibujas una función de densidad de probabilidad en R, dibujas una estimación de densidad tipo núcleo. El gráfico de densidad tipo núcleo es un enfoque no paramétrico que necesita seleccionar una ventana (bandwidth). Puedes establecer la ventana con el argumento bw de la función density.

En general, una ventana grande suavizará la curva de densidad, y una pequeña provocará un sobreajuste de la estimación de densidad tipo núcleo En el siguiente bloque de código encontrarás un ejemplo que describe este problema.

par(mfrow = c(1, 2))

# Ventana grande
plot(density(x, bw = 20), lwd = 2,
     col = "red", main = "Ventana demasiado grande")

# Ventana pequeña
plot(density(x, bw = 0.05), lwd = 2,
     col = "red", main = "Ventana demasiado pequeña")

De manera equivalente, puedes pasar argumentos de la función density a epdfPlot dentro de una lista como parámetro del argumento density.arg.list. En este caso, estamos pasando el argumento bw de la función de densidad.

# Alternativa equivalente con el paquete EnvStats
epdfPlot(x, epdf.col = "red", density.arg.list = list(bw = 0.05),
         main = "Ventana demasiado pequeña")

La literatura sobre la selección de la ventana para estimar funciones de densidad tipo núcleo es amplia. Sin embargo, hay tres enfoques principales comunmente utilizados para seleccionar el parámetro:

1. Por defecto, la función density utiliza el enfoque que se conoce como la “regla del pulgar”.
2. Usar el método plug-in, desarrollado por Sheather y Jones (1991).
3. Usar el enfoque de validación cruzada.

En el siguiente código se muestra cómo aplicar cada uno de los enfoques:

par(mfrow = c(1, 3))

# Regla del pulgar
plot(density(x), main = "Regla del pulgar",
     cex.lab = 1.5, cex.main = 1.75, lwd = 2)

# Validación cruzada insesgada
plot(density(x, bw = bw.ucv(x)), col = 2, # Igual a: bw = "UCV"
     main = "Validación cruzada", cex.lab = 1.5,
     cex.main = 1.75, lwd = 2)

# Plug-in
plot(density(x, bw = bw.SJ(x)), col = 4, # Igual a: bw = "SJ"
     main = "Método plug-in",
     cex.lab = 1.5, cex.main = 1.75, lwd = 2)

También puedes cambiar la función núcleo con el argumento kernel, que por defecto será gaussiana. Aunque no entraremos en más detalles, los núcleos disponibles son "gaussiano", "epanechnikov", "rectangular", "triangular", "biweight", "cosine" y "optcosine". La selección dependerá de los datos con los que estés trabajando.

Múltiples líneas de densidad en un gráfico

Con la función lines puedes dibujar múltiples curvas de densidad en R. Solo necesitas crear un gráfico de densidad en R y añadir las nuevas líneas que quieras.

par(mfrow = c(1, 1))

plot(dx, lwd = 2, col = "red",
    main = "Múltiples curvas", xlab = "")

# Nuevos datos
set.seed(2)
y <- rnorm(500) + 1
dy <- density(y)

# Nueva curva
lines(dy, col = "blue", lwd = 2)

Sin embargo, puede que te hayas dado cuenta de que la curva azul está cortada en el lado derecho. Para arreglar esto puedes establecer los argumentos xlim e ylim entre el mínimo y el máximo para cada eje de todas las densidades que vayas a dibujar.

plot(dx, lwd = 2, col = "red",
    main = "Multiples curvas con los límites de los ejes corregidos", xlab = "",
    xlim = c(min(dx$x, dy$x), c(max(dx$x, dy$x))),  # Mínimo y máximo limites eje X
    ylim = c(min(dx$y, dy$y), c(max(dx$y, dy$y))))  # Mínimo y máximo limites eje Y

lines(dy, col = "blue", lwd = 2)

Gráficos de comparación de densidades en R

Hay varias formas de comparar densidades. Un enfoque es utilizar la función densityPlot del paquete car. Esta función crea estimaciones de densidad no paramétricas condicionadas por un factor, si se especifica. Escribe ?densityPlot para obtener información adicional sobre la función y el método.

# Grupos
set.seed(1)
grupos <- factor(sample(c(1, 2), 100, replace = TRUE))

variable <- numeric(100)

# Grupo 1: media 3
variable[grupos == 1] <- rnorm(length(variable[grupos == 1]), 3)

# Grupo 2: media 0
variable[grupos == 2] <- rnorm(length(variable[grupos == 2]))

# Comparando densidades por grupo

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

densityPlot(variable, grupos)

Otra alternativa es utilizar la función sm.density.compare del paquete sm, que compara las densidades en un contraste de permutaciones.

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

sm.density.compare(variable, grupos)                
legend("topleft", levels(grupos), col = 2:4, lty = 1:2) 

Colorear el área bajo las curvas de densidad

Con R base puedes usar la función polygon para sombrear el área bajo las curvas de densidad. Si usas la función rgb en el argumento col en lugar de un color básico, puedes establecer la transparencia del área del diagrama de densidad con el argumento alpha, que va desde 0, que es transparencia total hasta 1, que dibuja el color opaco.

par(mfrow = c(1, 2))

#--------------------------------
# Sombrear el área bajo la curva
#--------------------------------

plot(dx, lwd = 2, main = "", xlab = "",
     col = "red", xlim = c(-4, 6), ylim = c(0, 0.5))
polygon(dx, col = "red")
polygon(dx$x, dx$y, col = "red") # Equivalente

set.seed(2)
y <- rnorm(500) + 2
dy <- density(y)

lines(dy, lwd = 2, col = "blue")
polygon(dy, col = "blue")


#--------------------------------------------------
# Sombrear el área bajo la curva con transparencia
#--------------------------------------------------

plot(dx, lwd = 2, main = "", xlab = "",
     col = "red", xlim = c(-4, 6), ylim = c(0, 0.5))
polygon(dx, col = rgb(1, 0, 0, alpha = 0.5))

lines(dy, lwd = 2, col = "blue")
polygon(dy, col = rgb(0, 0, 1, alpha = 0.5))

Si estás utilizando el paquete EnvStats, puedes configurar los colores con el argumento curve.fill.col de la función epdfPlot.

# Alternativa equivalente con el paquete EnvStats

library(EnvStats)

epdfPlot(x, # Vector con datos
        curve.fill = TRUE, # Colorear el área
        curve.fill.col = rgb(1, 0, 0, alpha = 0.5), # Color del área
        epdf.col = "red") # Color de la curva

epdfPlot(y, curve.fill = TRUE,
        curve.fill.col = rgb(0, 0, 1, alpha = 0.5),
        epdf.col = "blue",
        add = TRUE) # Añadir la densidad sobre el plot anterior

También puedes sombrear solo un área específica debajo de la curva. En el siguiente ejemplo, mostramos cómo colorear el área de la curva para valores de x mayores que 0.

par(mfrow = c(1, 1))

plot(dx, lwd = 2, main = "Densidad", col = "red")

polygon(c(dx$x[dx$x >= 0], 0), c(dx$y[dx$x >= 0], 0),
        col = rgb(1, 0, 0, alpha = 0.5), border = "red", main = "")

Gráfico de densidad con ggplot2

Puedes crear una gráfica de densidad con el paquete de R ggplot2 con las funciones ggplot y geom_density de la siguiente manera:

library(ggplot2)

df <- data.frame(x = x)

ggplot(df, aes(x = x)) + 
       geom_density(color = "red", # Color de la curva
                    fill = "red",  # Color del área sombreada
                    alpha = 0.5)   # Transparencia del color

Si quieres agregar más curvas, puedes establecer los límites del eje X con la función xlim y agregar una leyenda con scale_fill_discrete de la siguiente manera:

df <- data.frame(x = x, y = y)
df <- stack(df)

dx <- density(x)
dy <- density(y)

ggplot(df, aes(x = values, fill = ind)) + 
       geom_density(alpha = 0.5) + # Densidades con transparencia
       xlim(c(min(dx$x, dy$x), # Límites del eje X
            c(max(dx$x, dy$x)))) +
       scale_fill_discrete(name = "Título de la leyenda", # Cambiar el título de la leyenda
                           labels = c("A", "B")) # + # Cambiar las etiquetas de la leyenda
     # theme(legend.position = "none") # Eliminar leyenda


# Equivalente
ggplot(df, aes(x = values)) +
       geom_density(aes(group = ind, fill = ind), alpha = 0.5) + 
       xlim(c(min(dx$x, dy$x), c(max(dx$x, dy$x)))) +
       scale_fill_discrete(name = "Título de la leyenda", 
                           labels = c("A", "B"))