Dos gráficos en uno. R

Retomo mi olvidado blog para escribir algunas notas y para que no se me olvide como procedí cuando quería realizar un gráfico como el que presento a continuación.

Una de las cosas que mas me gusta de R es que primero uno piensa el gráfico que necesita, y luego lo hace, es decir, que los gráficos no están tan cerrados como en otros softwares estadísticos. En este caso la idea era realizar un gráfico de densidad de dos conjuntos de datos y añadirle un boxplot que coincidiera en el eje x. Esta idea es una modificación de la que en su día tuvo el Dr. Blas Benito para representar requerimientos ecológicos de especies vegetales y que se puede encontrar aquí.

Para comprenderlo mejor vamos a entrar en faena. Supongamos que tenemos un dataset con una variable binaria (presencia/ausencia) de una especie vegetal y muchas otras variables ambientales de dichos puntos.

Como dijimos nuestro propósito es representar gráficos de densidad de los datos de presencia y ausencia respecto a una variable ambiental concreta, y añadirle después un boxplot debajo, que coincida en el eje x. Para ello en primer lugar creamos dos subsets de datos (uno para los datos de presencia y otro para los datos de ausencia) (existen otras formas de hacerlo, pero yo he procedido de este modo):

# Separo los datos de presencia de los no presencia
datos0 <- datos[datos$presencia==0,]
datos1 <- datos[datos$presencia==1,] 

A continuación calculamos el gráfico de densidad para la variable ambiental, en este caso, la Temperatura mínima del invierno para el año 2000.

# Calculamos la Densidad para la variable temperatura mínima 
d.datos0 <- density(datos0$tmniob2000)
d.datos0 <- density(datos1$tmniob2000)

Seguidamente para asegurarnos de que todos los gráficos de todas las variables ambientales que queremos hacer guarden la misma proporción entre el plot de densidad y el boxplot, vamos a obtener los límites para los ejes x e y.

# Obtener los valores máximos de los ejes y el mínimo del eje x
ymax <- max(c(d.datos0$y, d.datos1$y))
xmin <- min(c(d.datos0$x, d.datos1$x))
xmax <- max(c(d.datos0$x, d.datos1$x))

El valor ymax representa el valor positivo máximo del eje y. Esta parte del gráfico (es decir, la parte positiva) queremos que siempre suponga el 70 % (por ejemplo) del gráfico combinado (boxplot y density). La parte negativa del eje y, será donde ubiquemos los boxplots y queremos que siempre ocupen un 30 % de la proporción del gráfico combinado. Para ello lo que vamos a hacer es crear una variable llamada ymin que sea el 30 % de la región del gráfico, teniendo en cuenta para ello el valor de ymax (El fijar estas proporciones es importante porque cuando realizamos este gráfico para diferentes variables ambientales, puede ser que los boxplots no tengan el mismo tamaño).

# Establecer un mínimo en el eje y, que sea el 30 % de la región del gráfico
ymin <- (ymax*(10/7))*0.3

Una vez que tenemos fijados dichos valores, vamos a empezar a realizar el gráfico. En primer lugar dibujamos la densidad de los datos de ausencias. Seguidamente añadimos la densidad de datos de presencia.

plot(d.datos0, ylim=c(-ymin,ymax), main="", ylab="", xlab="Temperatura mínima ivierno (ºC*10)")
polygon(d.datos1, col="#669700", border="#669700")
lines(d.datos0,xlim=c(0, max(datos0$tmxvob2000)), ylim=c(0,ymax), col="#336699", lwd=2) 

Tras hacerlo nos queda una cosa similar a esta

Si os fijáis nos queda un espacio vacío en la zona negativa del eje y, que siempre será de la misma proporción para todos los gráficos de todas las variables ambientales. Allí es donde vamos a realizar el gráfico de boxplots. Para ello vamos a utilizar una función llamada subplots(), del package TeachingDemos. Su funcionamiento es muy sencillo y como argumentos (ver ayuda de la función) establecemos qué gráfico vamos a incluir y en qué lugar. Para eso último utilizaremos las coordenadas de la región gráfica, haciendo uso del argumento 'usr' de la función par (). Este argumento es un vector de 4 elementos dando las coordenadas de los extremos de la región gráfica, en el siguiente modo c(x1, x2, y1, y2). Para el caso de las coordenadas del eje x, utilizamos par('usr')[1:2], y para la región negativa del eje y, lo que hacemos es dividirla en dos para que cuadren los dos boxplot. Todo ello lo hacemos así:

subplot(boxplot(datos1$tmxvob2000, horizontal=T, ylim=c(xmin,xmax), axes=F, col="#669700"),par('usr')[1:2], c((par('usr')[3] - 0)/2,0))
subplot(boxplot(datos0$tmxvob2000, horizontal=T, ylim=c(xmin,xmax), axes=F, col="#336699"),par('usr')[1:2], c(par('usr')[3],(par('usr')[3] - 0)/2))

Y conseguimos el gráfico combinado que nos proponíamos al inicio. Como se puede apreciar los boxplots coincíden en el eje x con los gráficos de densidad.

Pues nada más por ahora, espero os sirva.

Pd: Gracias al Dr. Benito por los datos.

Pd2: Para presentar el código he utilizado Pretty-R