El trabajo de cualquier analista de datos después de hacer una limpieza y consolidación de los mismos, es hacer una primer resumen de los mismos, bien sea gráfico o numérico. En este capítulo vamos a estudiar cómo podemos realizar la estadística descriptiva en R Commander. De este modo, se podrá tener una perspectiva general
Quiero leer...
Introducción a la Estadística Descriptiva en R
El objetivo principal de la estadística descriptiva es la de ordenar, resumir y clasificar los datos con el objetivo de tener una visión más precisa y global de las observaciones. De este modo, se pueden descubrir visualmente posibles relaciones entre los datos, ver sus puntos en común y sus diferencias, etc.
Pero otro punto importante de la estadística descriptiva es la de mostrar los datos de tal forma que permitan al investigador sugerir o proponer ciertas cuestiones para tratarlas posteriormente con mayor profundidad, así como determinar algunas suposiciones que son necesarias para poder realizar análisis estadísticos más complejos.
Matriz de datos
Lo primero que tenemos que tener para poder hacer estadística descriptiva son los datos. Y antes de empezar a ordenar, resumir y clasificar, tenemos que estructurar los datos en una matriz de datos (o una base de datos). Esta matriz se representa en modo de una tabla (Tabla 1), compuesta por columnas, que reflejan las variables del estudio y filas, que representan a los individuos sobre los que han sido medidas esas variables.
| Variable 1 | Variable 2 | … | Variable p | |
|---|---|---|---|---|
| Individuo 1 | x | x | … | x |
| Individuo 2 | x | x | … | x |
| Individuo 3 | x | x | … | x |
| … | … | … | … | … |
| Individuo n | x | x | … | x |
Veamos un ejemplo real de matriz de datos. Los datos que se exponen provienen de un artículo en el que medían volúmenes y algunas distancias particulares en los molares de Neandertales (Olejniczak et al. 2008). Como podemos observar en la Tabla 2 de ese estudio, las variables Enamel volume (mm3), Average enamel thickness (mm) y Relative enamel thickness (Volumen de esmalte, Grosor medio de esmalte y Grosor relativo de esmalte, respectivamente) se representan en columnas, mientras que los Neandertales (La Quina, Regourdou, Scladina…) en filas.
| Enamel volume (mm3) | Average enamel thickness (mm) | … | Relative enamel thickness | |
|---|---|---|---|---|
| La Quina | 222,11 | 1,15 | … | 14,67 |
| Regourdou | 171,16 | 0,93 | … | 13,14 |
| … | … | … | … | … |
| Scladina | 224,43 | 1,08 | … | 16,41 |
Representación gráfica de datos unidimensionales
La representación gráfica de los datos de una variable (unidimensionales) depende principalmente del tipo de variable que estemos usando. Para ello, recomendamos repasar los tipos de variables en el capítulo anterior, ya que cada tipo requiere un tipo de figura.
Gráficos de variables cualitativas
Veamos un ejemplo para poder representar gráficamente sus datos (Tabla 3). En un estudio sobre la expresión de las crestas del trigónido de molares de los homininos de la Sima de los Huesos de Atapuerca, se hizo el recuento sobre los cuatro tipos descritos de la variable cualitativa OES (Outer Enamel Surface, superficie external del esmalte) (Martínez de Pinillos et al. 2014). Los datos pueden descargarse desde la web del libro con el nombre molares.xlsx, y les llamamos en R Commander como Molares_SH. Este conjunto de datos tiene una columna llamada TipoOES.
| OES | Número | Frecuencia | Número acumulado | Frecuencia acumulada |
|---|---|---|---|---|
| A | 6 | 0.33 | 6 | 0.33 |
| B | 2 | 0.11 | 8 | 0.44 |
| C | 9 | 0.5 | 17 | 0.94 |
| D | 1 | 0.06 | 18 | 1 |
Se puede observar que se han descrito 4 tipos de OES (A, B, C, D) sobre un total de 18 molares de la población de la Sima de los Huesos.
Cálculo de frecuencias en R Commander
El cálculo de las frecuencias en R Commander es muy sencillo, aunque es necesario mencionar que en este caso, las frecuencias se muestran en tantos por 100, y no en base a 1. Observando la Figura 2 se comprenden cuáles son los pasos necesarios para calcularlo, ya que la ruta a seguir es Statistics - Summaries - Frequency distributions....
El resultado que nos devuelve R Commander es el siguiente:
.Table <- table(Molares_SH$TipoOES) .Table # counts for TipoOES
## ## A B C D ## 6 2 9 1
round(100*.Table/sum(.Table), 2) # percentages for TipoOES
## ## A B C D ## 33,33 11,11 50,00 5,56
Como puede observarse, primero aparece un recuento del número de individuos (en este caso molares) por categoría, y en segundo lugar ofrece la frecuencia relativa de los mismos en base 100.
Los datos de frecuencias de caracteres cualitativos se pueden representar gráficamente de dos modos: diagramas de barras y diagramas de sectores Para ver cómo se pueden hacer en R Commander, recomendamos ver la (Figura 3), siguiendo la ruta Graphs - Bar graph... / Pie chart..., respectivamente.
Diagramas de sectores en R Commander
Esta representación gráfica consiste en dividir un círculo en tantos sectores como categorías presentes haya en la variable cualitativa. En el caso del ejemplo, la variable cualitativa (TipoOES) tiene 4 categorías (A, B, C, D), por lo que el círculo presentará esas 4 categorías (Figura 10). La amplitud de cada sector depende del número o frecuencia de cada categoría. En R Commander se puede realizar siguiendo la ruta Graphs - Pie chart... (Figura 3 y Figura 4):
library(colorspace, pos=18)
with(Molares_SH, pie(table(TipoOES), labels=levels(TipoOES),
xlab="", ylab="", main="TipoOES", col=rainbow_hcl(4)))
Diagramas de barras en R Commander
Al igual que en el diagrama de sectores, un diagrama de barras consiste en construir tantos rectángulos como categorías presentes haya en la variable cualitativa. En el caso del ejemplo, la variable cualitativa (TipoOES) tiene 4 categorías (A, B, C, D), por lo que existirán 4 barras o rectángulos, cada uno por cada una de las cuatro categorías. La altura de cada barra depende del número o frecuencia de cada categoría. En R Commander se puede realizar siguiendo la ruta Graphs - Bar graph... (Figura 3 y Figura 5):
with(Molares_SH, Barplot(TipoOES, xlab="TipoOES", ylab="Frequency"))
Gráficos de variables cuantitativas
Los datos de variables cuantitativas pueden representarse dependiendo de si se encuentran agrupados en intervalos o no. Debido a que la representación de datos agrupados en intervalos es raro y exótico en estadística, no se abordará en el presente manual.
Sin embargo, es mucho más común y normal representar gráficamente los datos de variables cuantitativas que no están agrupadas en intervalos. Los gráficos principales son el diagrama de barras (si son pocos los valores distintos) o el histograma (si son muchos valores distintos).
Diagrama de caja y bigotes
Este diagrama de caja y bigotes se conoce también como diagrama de caja a secas o boxplot (de su procedencia anglosajona). Este diagrama está completamente relacionado con las medidas de tendencia central y con las medidas de dispersión.
Pongamos un ejemplo en el que se tomaron las medidas de anchura (diámetro mesiodistal, MD) y grosor (diámetro mesiobucal, BL) de los premolares inferiores de Neandertales (Lumley and Giacobini 2013). Nos vamos a centrar en concreto en los datos de la tabla 9 de dicho artículo. Esta tabla puede descargarse desde la web del libro con el nombre de Medidas MD y BL de P4 inf de neandertales – Lumley.txt. Cargamos el conjunto de datos en R Commander (Figura 6), y se le da el nombre de Lumley. Posteriormente se realiza el diagrama de caja y bigotes en R Commander siguiendo la ruta Graphs - Boxplot... (Figura 7).
En la Figura 8 aparece representado el diagrama de caja y bigotes del diámetro mesiobucal de los premolares de Neandertales. La única modificación que existe sobre el diagrama general que se obtuvo usando R Commander es la adición del nombre de los parámetros que están ahí representados (valores mínimo y máximo, primer y tercer cuartil, mediana y outlier). También pueden representarse sobre este gráfico la media (aritmética, geométrica, intercuartílica, etc) y la desviación estándar.
Este diagrama, a pesar de su aparente simpleza, presenta una altísima información resumida en sus líneas. De hecho, se hará especial incidencia a él y su relación con la distribución normal en el capítulo 3. El código que se ha ejecutado de modo predeterminado ha sido el siguiente:
Boxplot( ~ MD.mm., data=Lumley, id.method="y")
Este código lo que hace es indicar qué variable se va a representar (MD.mm.), los datos en los que se encuentra (data = Lumley) e identifica con números a los datos anómalos o outliers (id.method = "y"). Sin embargo, puede personalizarse un poquito más editando el código anterior (Figura 9) del siguiente modo:
Boxplot(~ MD.mm., # variable a representar
data = Lumley, # nombre del conjunto de datos
id.method = "n", # no muestra el número del outlier
col = "lightblue", # color de fondo de la caja
main = "Diagrama de caja y bigotes", # título de la figura
ylab = "Diámetro mesiodistal (mm)", # nombre del eje y
horizontal = FALSE, # representación en horizontal/vertical
notch = TRUE, # representación en modo notch
varwidth = FALSE) # control de la anchura de la barra
# Esto dibuja la cruz donde se sitúa la media
points(mean(Lumley$MD.mm.), # se calcula la media
col = "red", # color del punto
pch = 4, # tipo de símbolo
cex = 2) # tamaño
Como se observa, se ha cambiado el color de fondo (col="lightblue"), se ha dado un título a la figura y editado la leyenda del eje Y, indicando que se represente en modo notch. Finalmente se ha añadido una cruz con la media de la muestra. Se recomienda la lectura del Anexo 3 para ver un poco más en detalle la personalización de las figuras.
Histogramas
Un histograma es uno de los gráficos más comunes para representar datos cuantitativos. Realmente representa la frecuencia a intervalos concretos. Para poder interiorizarlo, se va a continuar con el ejemplo de este capítulo sobre los dientes de Neandertal. Se va a representar un histograma en R Commander para la variable MD.mm. como se indica en la Figura 10, siguiendo la ruta Graphs - Histogram....
El histograma para esta variable se puede observar en la Figura 11. Se observa que en el eje X (abscisas) está representado el intervalo de medidas de esa variable (entre 40 y 90 mm). Los valores que más se repiten (con una frecuencia de 11) son los comprendidos entre 70 y 75 mm.
La observación detallada de un histograma es vital para poder modelizar su distribución (por ejemplo una distribución normal), como se verá en detenimiento en el capítulo 3. El código que se ha ejecutado ha sido el siguiente:
with(Lumley, Hist(MD.mm., scale="frequency", breaks="Sturges", col="darkgray"))
Como siempre, se puede personalizar el código anterior modificándolo (Figura 12):
with(Lumley, Hist(MD.mm., scale="density", breaks=20, col="darkred"))
La observación detallada de un histograma es vital para poder modelizar su distribución (por ejemplo una distribución normal), como se verá en detenimiento en el capítulo 3. El código que se ha ejecutado ha sido el siguiente:
with(Lumley, Hist(MD.mm., scale="frequency", breaks="Sturges", col="darkgray"))
Como siempre, se puede personalizar el código anterior modificándolo (Figura 12):
with(Lumley, Hist(MD.mm., scale="density", breaks=20, col="darkred"))
En este código se ha representado la densidad (también se puede cambiar en la pestaña Opciones de la interfaz gráfica), el color y el número de barras (breaks = 20). Respecto a este último argumento, de modo predeterminado aparece Sturges, pero podemos indicarle manualmente cuántas barras queremos (como en el ejemplo, breaks=20). También puede utilizarse otra aproximación, como la de Freedman-Diaconis, teniendo que ponerlo de la siguiente forma (breaks = FD).
Curvas de densidad
Las curvas de densidad representan cómo se distribuyen los datos a lo largo de toda la variable. Es similar a un histograma, pero en vez de representarlo mediante barras, es solo una línea que muestra cómo se dispersan y se concentran los datos. Para poder comprenderlo, se continúa con el ejemplo de este capítulo sobre los dientes de Neandertal. Una curva de densidad se ejecuta en R Commander, para la variable MD.mm., como se indica en la Figura 13, siguiendo la ruta Graphs - Density estimate....
El resultado se muestra en la Figura 14. Como se puede observar, aparte de la curva de densidad (que tiene una especie de forma de montaña), se observan unas pequeñas líneas verticales situadas justo por encima del eje X. Estas líneas representan las observaciones (en nuestro caso, las medidas reales de los dientes de Neandertal). Cuanto más juntas estén, más alta es la curva. Este es el código que se ha ejecutado automáticamente para representar la curva de densidad:
densityPlot( ~ MD.mm., data=Lumley, bw="SJ", adjust=1, kernel="gaussian")
Representación gráfica de datos bidimensionales
En este caso se va a representar uno de los diagramas más conocidos en la estadística: el diagrama de dispersión (o scatterplot en inglés). En él, dos variables se enfrentan, una representada en el eje X (eje abscisas) y otro en el eje Y (eje de ordenadas).
Siguiendo con el ejemplo anterior de los diámetros de los premolares de Neandertal, se va a realizar un diagrama de dispersión usando R Commander (Figura 15), mostrando en el X el diámetro bucolingual (BL.mm.) y en el eje Y el diámetro mesiodistal (MD.mm.). La ruta a seguir en R Commander es Graphs - Scatterplot....
Observando la Figura 16 se puede adivinar a simple vista, que cuando aumenta el diámetro bucolingual, tiende a aumentar también el mesiodistal. Aunque en este capítulo no se profundizará en cómo es esa relación, en el capítulo 5 se estudiará en profundidad la regresión y correlación de datos bivariantes.
scatterplot(MD.mm.~BL.mm., reg.line=FALSE, smooth=FALSE, spread=FALSE, boxplots='xy', span=0.5, ellipse=FALSE, levels=c(.5, .9), data=Lumley)
Diagrama de dispersión dividido por grupos
Pero vamos a complicar un poco más el proceso. Vamos a añadir una columna a la derecha haciendo que aparezca como se ve en la Tabla 4. Los datos os los podéis descargar desde la web del libro, con el nombre de Lumley con grupos en columna especie.txt. En este conjunto de datos el símbolo decimal es un punto. Cargamos el conjunto de datos nuevamente en R Commander y le damos el nombre de Lumley_sp (o cualquier otro, a elección).
| P4.inferior | MD.mm. | BL.mm. | especie |
|---|---|---|---|
| Fate 2 g | 7,3 | 8,1 | nean.fate |
| Fossellone 3g | 7,4 | 8,5 | nean.fossellone |
| Zafarraya 2 d | 7,0 | 9,5 | nean.zafarraya |
| Zafarraya 2 g | 7,1 | 9,0 | nean.zafarraya |
| Genay d | 7,5 | 9,8 | nean.genay |
| Genay g | 7,4 | 10,1 | nean.genay |
| Le Moustier | 8,3 | 10,2 | nean.moustier |
| Arcy II | 7,8 | 9,5 | nean.arcy |
| Regourdou | 6,1 | 8,6 | nean.regourdou |
| La Quina 5 d | 6,6 | 10,5 | nean.quina |
| La Quina 5 g | 7,9 | 9,2 | nean.quina |
| Krapina min. | 7,2 | 9,2 | nean.krapina |
| Krapina max. | 8,6 | 11,3 | nean.krapina |
| Hortus II d | 8,0 | 8,4 | nean.hortus |
| Hortus II g | 7,8 | 8,4 | nean.hortus |
| Hortus V g | 7,0 | 8,2 | nean.hortus |
| Mauer d | 7,5 | 9,2 | nean.mauer |
| Arago 2 g | 7,3 | 9,3 | nean.arago |
| Arago 13 d | 8,5 | 11,5 | nean.arago |
| Arago 28 g | 8,8 | 10,2 | nean.arago |
| Atapuerca 2 d | 7,5 | 8,0 | nean.atapuerca |
| Atapuerca 2 g | 6,5 | 8,2 | nean.atapuerca |
| Atapuerca 3 d | 7,4 | 10,2 | nean.atapuerca |
| Atapuerca 6 g | 8,2 | 10,5 | nean.atapuerca |
| Bañolas d | 4,9 | 7,8 | nean.bañolas |
| Bañolas g | 6,0 | 8,2 | nean.bañolas |
| Hommes actuels min. | 6,5 | 7,0 | sapiens |
| Hommes actuels moy. | 7,1 | 8,0 | sapiens |
| Hommes actuels max. | 8,0 | 9,0 | sapiens |
Ahora se representa un diagrama de dispersión en el que se diferencien los puntos según el grupo que aparezcan en la variable cualitativa especie. Para ello, siguen los pasos de la Figura 15, aunque con la salvedad de hacer presión sobre el botón Plot by groups..., tal y como se muestra en la Figura 17.
El resultado de este diagrama de dispersión se puede observar en la Figura 18. Manteniendo las mismas opciones que en la pestaña Options de la Figura 15, se puede conseguir fácilmente este diagrama de dispersión.
scatterplot(MD.mm.~BL.mm. | especie, reg.line=FALSE, smooth=FALSE, spread=FALSE, boxplots=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), by.groups=TRUE, data=Lumley_sp)
Una funcionalidad muy interesante es la de representar las elipses de equiprobabilidad sobre el diagrama de dispersión, bien para toda la muestra en conjunto, bien por cada grupo. En la Figura 15 se observa una opción para marcar que dice Plot concentration ellipse(s). De modo automático te marca que se hagan dos elipses, una con una probabilidad de 0,5 y la otra 0,9. Esto es, obviamente, personalizable. En la Figura 19 se muestra el resultado de estas dos elipses, junto con su centroide (punto verde). También se pueden dibujar teniendo en cuenta cada uno de los grupos, tal y como se observa siguiendo los pasos de la Figura 17.
scatterplot(MD.mm.~BL.mm., reg.line=FALSE, smooth=FALSE, spread=FALSE, boxplots=FALSE, span=0.5, ellipse=TRUE, levels=c(.5, .9), data=Lumley_sp)
Matrices de diagramas de dispersión
Cuando se tiene más de dos variables numéricas y se quiere intentar interpretar cómo están relacionadas entre ellas mediante un diagrama de dispersión, lo más sensato es representarlas enfrentándose una a una. Esto se consigue con una matriz de diagramas de dispersión. Se va a emplear los datos llamados litica_silex.xlsx y que pueden descargarse desde la web del libro. Estos datos representan un estudio sobre la caracterización del sílex sobre el que se han tomado algunas variables (IRA, IRC y Peso) que ayudan a identificar su procedencia y contextualizar movimientos de aprovisionamiento de esta materia prima durante el Pleistoceno y el Holoceno. Se unoirta el archivo excel siguiendo la siguiente ruta en R Commander Data - Import data - from Excel file....
En R Commander seguimos la siguiente ruta, indicando que nos distinga los grupos (Figura 20): Graphs - Scatterplot matrix....
En la Figura 21 se observa el resultado. En la diagonal aparecen las curvas de densidad, pero pueden ser cambiadas por otro tipo de figuras, como los diagramas de cajas y bigotes.
scatterplotMatrix(~IRA+IRC+Peso | Grupo.Laplace, reg.line=FALSE, smooth=FALSE, spread=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), id.n=0, diagonal= 'density', by.groups=TRUE, data=litica)
Áreas de densidad con líneas
Para representar gráficamente áreas de densidad se necesita instalar y cargar el paquete ggplot2 con las funciones install.packages("ggplot2") y library(ggplot2). Su función principal es la representación gráfica y contiene prácticamente infinitas opciones de visualización.
En una publicación reciente se han tomado las medidas de longitud de dos dedos (ANULAR e ÍNDICE) y de la mano (LARGO) en 46 mujeres y en 31 hombres, con el objetivo de que sirvan de referencia para compararlas con las manos paleolíticas de la Cueva de El Castillo (Cantabria, España) (Rabazo-Rodríguez et al. 2017).
Los datos del ejemplo pueden descargarse desde la web del libro, con el nombre medidas_manos.txt. En esta ocasión, se va a llamar al conjunto de datos como manos. Se emplea la variable LARGO como nuestra variable en el eje X y la variable ANULAR como la homóloga del eje Y. Siguiendo los pasos descritos en la sección Diagrama de dispersión dividido por grupos, indicando que la variable SEXO sea la variable categórica, se obtiene la Figura 22.
scatterplot(ANULAR~LARGO | SEXO, reg.line=FALSE, smooth=FALSE, spread=FALSE, boxplots=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), by.groups=TRUE, data=manos)
Como se puede observar, es un diagrama de dispersión sencillo, sin excesos estilísticos. Es informativo, pero su visualización es común, pudiendo existir ciertos aspectos que no sean fácilmente identificables con su escrutinio visual, como por ejemplo, ver cómo cambia la densidad de sus puntos, desde el área con mayor concentración hacia la zona con menor densidad.
En la Figura 23 se muestran dos gráficos, el primero representando el área de densidad de todos los datos en conjunto (A), y el segundo discriminando por SEXO (B).
library(ggplot2) # cargamos el paquete # Código para representar las curvas de densidad sin tener en cuenta SEXO ggplot(manos, aes (x = LARGO, y = ANULAR)) + stat_density2d () # Código para representar las curvas de densidad teniendo en cuenta SEXO ggplot(manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) + stat_density2d ()
Los dos códigos anteriores están divididos en dos partes, aunque conectados utilizando el símbolo más (+) como separador:
ggplot(). En esta función hay que añadir la información de los datos, indicando cómo se llama el conjunto de datos (manos), la variable X (LARGO), la variable Y (ANULAR) y la variable que discrimina por grupos (colour = SEXO).stat_density2d(). Esta función representa gráficamente las curvas de densidad. De modo predeterminado ejecuta el argumentogeom = "contour". Esto quiere decir que ejecutarstat_density2d()ystat_density2d(geom = "contour")va a dar el mismo resultado.
Personalizar la función stat_density2d() mediante la inclusion de nuevos argumentos permite mejorar la visualización de la figura. Además, como puede verse en la Figura 23, algunas curvas de densidad, más específicamente las externas, están cortadas en los bordes del gráfico. Este hecho hace necesario añadir nuevas funciones a que controlen los límites numéricos de los ejes X e Y, tal y como se ven en los dos códigos siguientes y su representación gráfica en la Figura 24.
ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) + stat_density2d (lty = 2, lwd = 0.5) + xlim(140, 220) + ylim(62, 92) ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) + stat_density2d (lty = 2, lwd = 0.5, bins = 20) + xlim(140, 220) + ylim(62, 92)
Conociendo los argumentos lty y lwd, explicados en el anexo 3, se pueden personalizar las líneas de las curvas de densidad dentro de stat_density2s(). Por otro lado, el argumento bins = 20 precisa el número de curvas de densidad por cada grupo, siendo en este caso 20. Las dos nuevas funciones que se han incluido para personalizar los intervalos de los ejes X e Y son las siguientes:
xlim(). Se controla el intervalo del eje X.ylim(). Se controla el intervalo del eje Y.
Se pueden continuar personalizando los códigos, como se ven en los que se encuentran a continuación y están representados gráficamente en la Figura 25.
ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) +
stat_density2d (aes (fill = SEXO, colour = SEXO, alpha = 0.9),
geom = "polygon", lty = 3, lwd = 1) +
xlim(140, 220) +
ylim(62, 92) +
scale_alpha(range = c(0, 1/2), guide = "none")
ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) +
stat_density2d (aes (fill = SEXO, colour = SEXO, alpha = 0.9),
geom = "polygon", lty = 3, lwd = 1) +
xlim(140, 220) +
ylim(62, 92) +
scale_alpha(range = c(0, 1/2), guide = "none") +
geom_point(aes(fill = SEXO), colour = "black", pch = 21, size = 3, alpha = 0.5)
Aquí se empieza complicar un poquito más el código, pero vamos a explicarlo para que conozcáis perfectamente qué es lo que se ha ejecutado para generar las figuras. Analizaremos igual las otras secciones incluidas:
stat_density2d(). Se ha incluidogeom = "polygon", que pinta el área de las curvas de colores. Para que las pinte de los mismos colores que las curvas, tenemos que añadir lo especificado enaes(). Dentro de esta última función hay que especificar cuál es la variable categórica (SEXO), así como la transparencia (alpha).scale_alpha(). Esta función relativiza los colores. Me explico, si no lo incluyéramos, aparecería un manchurrón azul o rojo, dependiendo del grupo, en la zona central, y apenas se vería transición entre la zona más densa a la menos densa. Por lo que su función es garantizar que visualmente exista ese degradado.geom_point(). Esta sección añade los puntos de nuestras observaciones. Igual que en el primer caso, hay que especificar cuál es la variable categórica (SEXO) para que tengan los mismos colores. El resto de argumentos están especificados en el anexo 3.
Echamos un vistazo a los 6 códigos siguientes, representados gráficamente en la Figura 26:
# Código 1A ggplot (manos, aes (x = LARGO, y = ANULAR)) + stat_density2d (geom = "point", aes(size = ..density..), n = 20, contour = FALSE) # Código 1B ggplot (manos, aes (x = LARGO, y = ANULAR, colour=SEXO)) + stat_density2d (geom = "point", aes(size = ..density..), n = 20, contour = FALSE) # Código 2A ggplot (manos, aes (x = LARGO, y = ANULAR)) + stat_density2d (geom = "raster", aes(fill = ..density..), contour = FALSE) # Código 2B ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) + stat_density2d (geom = "raster", aes(fill = ..density..), contour = FALSE) # Código 3A ggplot (manos, aes (x = LARGO, y = ANULAR)) + stat_density2d (aes(fill = ..level..), geom = "polygon") # Código 3B ggplot (manos, aes (x = LARGO, y = ANULAR, colour = SEXO)) + stat_density2d (aes(fill = ..level..), geom = "polygon")
En la Figura 26 se han representado algunas formas complementarias de mostrar las densidades de nuestros datos en un diagrama de dispersión.
En los códigos 1A y 1B se muestran cómo se vería representando por puntos de densidad. Cuanto más grande es el punto, más concentración de datos existe. Destaca la figura 1b por representar gráficamente los puntos que contienen datos de los dos grupos. Los puntos más pequeños de las esquinas superior izquierda e inferior derecha indican ausencia de datos.
La representación mostrada en los códigos 2A y 2B son también interesantes. Cuanto más claro sea el color azul, más densidad. Aunque en la figura 2a es fácil interpretarlo, no es así en la 2b, donde no se permite diferenciar los dos grupos con claridad.
El último caso, que agrupa los códigos 3A y 3B, se representa el relleno simplificado de la Figura 30. Es más, en la 3b, como algunas líneas de la curva están cortadas, el relleno de color genera unas anomalías que se solucionarían sin problemas introduciendo los límites de los ejes X e Y, como se observan en la Figura 24 y Figura 25.
Áreas de densidad con hexágonos
Para finalizar esta sección dedicada a la representación gráfica de datos bivariantes, vamos a aprender a representar gráficamente un diagrama dispersión mostrando las densidades de nuestros datos mediante hexágonos de color. Este gráfico requiere los paquetes ggplot2 y hexbin2. Como el primero ya lo tenemos instalado y cargado (página 55), instalamos y cargamos el segundo:
install.packages("hexbin")
library(hexbin)
A continuación ejecutamos los siguientes códigos, aplicados al mismo ejemplo de las longitudes de los dedos y manos (Figura 27):
# Código A ggplot (manos, aes (x = LARGO, y = ANULAR)) + stat_binhex (bins=5, aes (alpha = ..count..)) # Código B ggplot (manos, aes (x = LARGO, y = ANULAR, fill = SEXO)) + stat_binhex (bins=5, aes (alpha = ..count..)) # Código C ggplot (manos, aes (x = LARGO, y = ANULAR, fill = SEXO)) + stat_binhex (bins=10, aes (alpha = ..count..))
Los códigos anteriores se corresponden con el mismo orden, de arriba a abajo, representado en la Figura 27. En este caso, observamos que la sección personalizada de los hexágonos se llama stat_binhex(). Mediante el atributo bins, podremos aumentar o disminuir el número de hexágonos.
Si queremos tener separados los grupos de la figura, apareciendo en nuestro ejemplo los hombres a un lado y las mujeres a otro, tenemos que personalizar los códigos anteriores añadiendo a cada uno esta sección, facet_grid (. ~ SEXO), tal y como se observa a continuación (Figura 28):
# Código A ggplot (manos, aes (x = LARGO, y = ANULAR)) + stat_binhex (bins=5, aes (alpha = ..count..)) + facet_grid (. ~ SEXO) # Código B ggplot (manos, aes (x = LARGO, y = ANULAR, fill = SEXO)) + stat_binhex (bins=5, aes (alpha = ..count..)) + facet_grid (. ~ SEXO) # Código C ggplot (manos, aes (x = LARGO, y = ANULAR, fill = SEXO)) + stat_binhex (bins=10, aes (alpha = ..count..)) + facet_grid (. ~ SEXO)
Representación gráfica de datos tridimensionales
Para representar gráficamente datos tridimensionales en R es necesario instalar unos paquetes extras. En esta sección vamos a aprender a representarlos usando dos paquetes: scatterplot3d y rgl. Y como ejemplo, vamos a representar en tres dimensiones diagramas de dispersión compuestos por tres variables. La ventaja del segundo respecto al primero es que es interactivo, es decir, puedes desplazar con el ratón el cubo generado para buscar la mejor visualización de los datos, mientras que en el primero tienes que indicarle los grados desde donde ver tus datos. Independientemente, ambos son una buena herramienta para ver en 3D la distribución de nuestros datos. Vamos a utilizar el conjunto de datos manos que presenta tres variables cuantitativas: INDICE, ANULAR y LARGO.
Diagramas de dispersión en 3D con scatterplot3d
scatterplot3d es un paquete relativamente sencillo de utilizar. Pero antes de ponernos manos a la obra con él tenemos que instalarlo y cargarlo en R Commander. Para ello ejecutamos los dos códigos siguientes:
install.packages("scatterplot3d")
library(scatterplot3d)
Con el conjunto de datos manos cargado, vamos a ejecutar la función scatterplot3d()(Figura 29):
scatterplot3d(manos$INDICE, manos$ANULAR, manos$LARGO)
En esta función hay que definir primero cuáles son nuestras variables X, Y y Z, respectivamente. Por lo tanto, hemos añadido cuáles son las variables X (manos$INDICE), Y (manos$ANULAR) y Z (manos$LARGO). Recordamos que el símbolo de dólar ($) lo que hace es unir el nombre del conjunto de datos (manos) con el nombre de la columna de una variable concreta (INDICE, ANULAR, LARGO).
Sin embargo, podemos personalizar la figura anterior. Atención al siguiente código a ejecutar (Figura 30):
# Parte 1. Como tenemos 2 en sexo, definimos qué símbolo representa cada uno.
shapes=c(14,15)
shapes=shapes[as.numeric(manos$SEXO)]
# Parte 2. Como tenemos 2 en sexo, definimos qué color representa cada uno.
colors=c("red", "blue")
colors=colors[as.numeric(manos$SEXO)]
# Parte 3. Incluimos más argumentos en la función.
scatterplot3d(manos$INDICE, manos$ANULAR, manos$LARGO,
pch = shapes,
type = "h",
color = colors,
grid = TRUE,
box = TRUE,
angle = 130,
main = "Diagrama de dispersión en 3D",
xlab = "X = Índice",
ylab = "Y = Anular",
zlab = "Z = Largo")
Vamos a intentar comprender el código anterior. Como véis, lo he dividido en 3 partes. Aun así, os comento un elemento clave en programación para que os vayáis familiarizando con él. Todo el código que está escrito después del símbolo almohadilla (#) es código que no se va a ejecutar. Por eso veréis en muchos sitios (y en este libro también) que se utiliza para hacer comentarios sobre el código, principalmente para facilitar su comprensión y saber qué va a hacer. Volvamos a las 3 partes del código:
- Parte 1: está escrita para personalizar los símbolos de cada grupo. Como tenemos 2 grupos en nuestros datos, tenemos que especificar 2 tipos de símbolos, uno por grupo. En este caso, se hace con números (14 y 15). En el anexo 3 podréis ver la equivalencia de los símbolos son sus respectivos números (
pch). - Parte 2. De modo análogo a la parte 1, aquí se especifican los colores de cada grupo. Igualmente, en el anexo 3 podréis encontrar información sobre los colores para que lo personalicéis con vuestras preferencias.
- Parte 3. Una vez especificados los puntos anteriores, vamos a ejecutar la función propiamente dicha de
scatterplot3d(). Pero vamos a modificar la función añadiendo más atributos:- Como ya se ha explicado, hay que definir inicialmente las variables X, Y, Z (
manos$INDICE,manos$ANULAR,manos$LARGO). pchycolorhacen referencia a los objetos creados en la Parte 1 y Parte 2, respectivamente.typese refiere a si se quieren representar solo los puntos (p), líneas que unen los puntos (l) o la línea que une cada punto con el plano inferior del 3D (h).gridyboxpueden adquirir los valoresTRUE / FALSE, y representan la malla cuadrada del plano inferior o el cubo 3D donde se ubican los puntos, respectivamente.main,xlab,ylab,zlab, se refieren a los nombres del título de la figura, y de los ejes X, Y, Z respectivamente.anglepermite cambiar el ángulo de la figura 3D.
- Como ya se ha explicado, hay que definir inicialmente las variables X, Y, Z (
Sin embargo, y esto ya es una opinión personal, me resulta complicado poder interpretar los datos utilizando esta figura, ya que no te permite tener libertad de rotar el 3D para poder facilitar la comprensión de los datos. Pero el que sí te lo permite es el paquete que vamos a aprender a continuación (rgl).
Diagramas de dispersión en 3D con el paquete rgl
El paquete rgl es más completo y complejo que scatterplot3d. Además te permite generar figuras interactivas que puedes mover interaccionando con el ratón. De entre las muchas opciones que posee a la hora de representar gráficamente se ha seleccionado un tipo de visualización que permite representar las observaciones en 3D, diferenciando por grupos y marcando las elipses de equiprobabilidad cada uno de ellos. Se puede acceder a esta opción en R Commander siguiendo la ruta Graphs - 3D graph - 3D scatterplot..., donde nos aparece la ventana de la Figura 31. En esta ventana se seleccionan dos variables explicativas en la izquierda y una variable respuesta en la derecha. La ejecución de las funciones de esta ventana automáticamente cargan el paquete rgl.
Con las opciones predeterminadas y habiendo seleccionado ANULAR e INDICE como variables independientes y LARGO como variable dependiente, se obtiene la Figura 32.
scatter3d(LARGO~ANULAR+INDICE, data=manos, surface=FALSE, residuals=TRUE, bg="white", axis.scales=TRUE, grid=TRUE, ellipsoid=FALSE)
Si además ahora se desea dividiar por grupos y mostrar las elipses de equiprobabilidad, en la ventana de la Figura 31 se selecciona Plot by groups..., donde se marca la variable SEXO como variable categórica, y la pestaña de Options se marca la opción de Plot 50% concentration ellipsoid. Con estas opciones establecidas se obtiene la Figura 33
scatter3d(LARGO~ANULAR+INDICE|SEXO, data=manos, surface=FALSE, residuals=TRUE, parallel=FALSE, bg="white", axis.scales=TRUE, grid=TRUE, ellipsoid=TRUE)
Medidas de tendencia central
Al margen de las figuras y gráficos estudiados anteriormente, también es posible representar numéricamente información útil mediante números. Estos números resumen un conjunto de datos. En conjunto se conocen como promedios, medidas de posición o medidas de tendencia central. Entre otros, algunas medidas de tendencia central son la media aritmética, la mediana, la moda y los cuantiles.
Continuando con el ejemplo de los dientes de Neandertales (Lumley and Giacobini 2013), se han ordenado de menor a mayor todos los valores de los diámetros mesiodistales de sus dientes (Figura 34).
La media aritmética se calcula sumando todos los valores de una muestra y dividiéndolo por el número total de la muestra, mientras que la moda es el valor que más se repite de los datos. Si varios valores presentan igual frecuencia, la moda se comparte entre ellos. Es decir, habrá tantas modas como datos con la frecuencia más alta e igual.
Los cuantiles son una medida de posición. Se expresa con el cuantil y un valor entre 0 y 1. Estos valores se corresponden con una posición concreta de las observaciones de nuestra variable cuando los datos son ordenados de mayor a menor. Si ordenamos todas las observaciones de menor a mayor, la mediana es el valor que se encuentra justo en el medio. Este valor se corresponde por tanto con el cuantil 0.50. Si el número de observaciones es par, la mediana se calcula haciendo la media de los dos valores centrales. Por ejemplo, si tenemos los valores siguientes 2.2, 3.1, 2.6, 3.2 y 2.9 (en total, 5 datos), lo primero que hacemos es ordenadors de menor a mayor: 2.2, 2.6, 2.9, 3.1 y 3.2. El cuantil 0.5 (o mediana) es el valor que se encuentra justo en la mitad, en este caso siendo 2.9.
Los cuartiles son un tipo especial de cuantil en los que se representa el 0.25, 0.50 y 0.75, es decir, la muestra se divide en cuartos (de ahí el nombre de cuartiles), y se da el valor de nuestros datos que se corresponderían con la posición del cuantil 0.25, 0.50 y 0.75 (Figura 34). El rango intercuartílico es el que se produce de la diferencia entre el valor del cuartil 0.75 menos el cuartil 0.25. Por lo tanto, todos los cuartiles son cuantiles, pero no todos los cuantiles son cuartiles. Se podría entender de todo esto, además, que el cuantil 0.00 se corresponde con el valor mínimo, el cuantil 0.50 con la mediana y el cuantil 1.00 con el valor máximo.
En este ejemplo se ha medido el diámetro mesiodistal sobre 29 premolares de Neandertal. Lo que importa para comprender los cuantiles es ese número, el 29. Ahora, lo que hacemos es, siguiendo lo descrito previamente, ordenar todos los datos de menor a mayor. En R Commander se realiza siguiendo la ruta Statistics - Summaries - Numerical summaries... (Figura 35).
Se despliega a continuación una ventana con dos pestañas (Data | Statistics). En la primera de ellas, se seleccionan la/s variable/s de las que se quiere obtener las medidas de tendencia central, mientras que en la segunda se seleccionan qué medidas exactamente se desea obtener, como la media, desviación estándar, rango intercuartílico y cuantiles (0.00, 0.25, 0.50, 0.75, 1.00).
Se presiona sobre OK y se obtiene el siguiente resultado para la variable MD.mm.:
numSummary(Lumley[,"MD.mm.", drop=FALSE], statistics=c("mean", "sd", "IQR",
"quantiles"), quantiles=c(0,.25,.5,.75,1))
## mean sd IQR 0% 25% 50% 75% 100% n ## 73,51724 8,458388 9 49 70 74 79 88 29
Una alternativa más rápida, pero con algunas opciones menos, es la ofrecida por R Commander y que puede ejecutarse al seguir la ruta Statistics - Summaries - Active data set.... En este caso, aparecen automáticamente los cuartiles y la media de todas las variables de nuestros datos (en este caso, MD.mm. y BL.mm.), como se ve a continuación:
summary(Lumley)
## P4 inferior MD.mm. BL.mm. ## Length:29 Min. :49,00 Min. : 70,00 ## Class :character 1st Qu.:70,00 1st Qu.: 82,00 ## Mode :character Median :74,00 Median : 92,00 ## Mean :73,52 Mean : 91,59 ## 3rd Qu.:79,00 3rd Qu.:101,00 ## Max. :88,00 Max. :115,00
Medidas de dispersión
Se ha estudiado en la sección anterior que las medidas de tendencia central te ofrecen un valor numérico que resume los datos, como la media aritmética. Ahora bien, ¿están muy concentrados los datos en torno a esa media? O por el contrario… ¿están muy separados o dispersos?
Para responder a estas preguntas se recurre a las medidas de dispersión. El objetivo de las medidas de dispersión es, por tanto, el de estudiar lo concentrados que están los datos en torno a un promedio. Cuanto más bajo sea el valor de dispersión, más concentrados estarán en torno a la media. Entre otros, algunas medidas de dispersión son la varianza y la desviación típica.
Conviene comentar que la varianza y la desviación típica (también conocida como desviación estándar) pueden obtenerse de la población completa de la que queramos obtener información, o bien de la muestra que tengamos de esa población. Las fórmulas para su cálculo son diferentes dependiendo de si es la población o una muestra. Para ser técnicamente correctos en el empleo de los conceptos, el cálculo de ambas medidas de dispersión sobre las muestras se denominan cuasivarianzas y cuasidesviaciones típicas, respectivamente.
Como prácticamente siempre se trabaja con muestras (no sólo en arqueología, sino también en el resto de disciplinas científicas), realmente se utilizarán las fórmulas de las cuasivarianzas y cuasidesviaciones. Pero es tan común el empleo de muestras en los estudios científicos que prácticamente se ha desechado el uso del prefijo cuasi- al referirse a esas medidas de dispersión. Incluso esto sucede en R, que se refiere a las cuasivarianzas y cuasidesviaciones como varianzas y desviaciones respectivamente.
Existe una relación matemática entre la desviación típica y la varianza. La desviación típica es la raíz cuadrada de la varianza, o lo que es lo mismo, la varianza es la desviación típica elevada al cuadrado.
Se ha medido en un estudio de la extensión en hectáreas de 25 asentamientos prehistóricos del Uruk tardío, en la antigua Mesopotamia (Johnson 1973). Este conjunto de datos puede descargarse desde la web del libro con el nombre de Datos Uruk.txt. Se le da el nombre de Uruk.
Una vez cargados los datos en R Commander, y habiendo seguido la ruta Statistics - Summaries - Numerical summaries... (Figura 36), se puede obtener la desviación típica (sd de standard deviation). Como se observa, el valor de la desviación estándar es 18.73 hectáreas.
numSummary(Uruk[,"Superficie", drop=FALSE], statistics=c("mean", "sd",
"IQR", "quantiles"), quantiles=c(0,.25,.5,.75,1))
## mean sd IQR 0% 25% 50% 75% 100% n ## 52,124 18,72563 28 30 37 45 65 90,5 25
El cálculo de la varianza hay que hacerlo manualmente mediante la ejecución de la función var(). En este caso, la función var() se refiere a la varianza, y entre paréntesis se escribe el nombre de los datos activos + el símbolo de dólar + el nombre de la variable de la que queremos calcular la varianza. Se ejecuta el código y se observa que la varianza es 350.65 hectáreas.
var(Uruk$Superficie)
## [1] 350,6494
Referencias
Olejniczak, Anthony J, Tanya M Smith, Robin N M Feeney, Roberto Macchiarelli, Arnaud Mazurier, Luca Bondioli, Antonio Rosas, et al. 2008. “Dental Tissue Proportions and Enamel Thickness in Neandertal and Modern Human Molars.” Journal of Human Evolution 55 (1): 12–23. http://www.sciencedirect.com/science/article/B6WJS-4S094W2-1/2/967306ae20e21a9495032d2724941a35.
Martínez de Pinillos, Marina, María Martinón-Torres, Matthew M. Skinner, Juan Luis Arsuaga, Ana Gracia-Téllez, Ignacio Martínez, Laura Martín-Francés, and José María Bermúdez de Castro. 2014. “Trigonid Crests Expression in Atapuerca-Sima de Los Huesos Lower Molars: Internal and External Morphological Expression and Evolutionary Inferences.” Comptes Rendus Palevol 13 (3): 205–21. https://doi.org/10.1016/j.crpv.2013.10.008.
Lumley, Marie-Antoinette de, and Giacomo Giacobini. 2013. “Les Néandertaliens de La Caverna Delle Fate (Finale Ligure, Italie). II Les Dents.” L’Anthropologie 117 (3): 305–44. https://doi.org/10.1016/j.anthro.2013.05.002.
Rabazo-Rodríguez, Ana María, Mario Modesto-Mata, Lucía Bermejo, and Marcos García-Díez. 2017. “New Data on the Sexual Dimorphism of the Hand Stencils in El Castillo Cave (Cantabria, Spain).” Journal of Archaeological Science: Reports 14 (June): 374–81. https://doi.org/10.1016/j.jasrep.2017.06.022.
Johnson, Gregory Alan. 1973. Local Exchange and Early State Development in Southwestern Iran. 1st ed. Michigan: Museum of Anthropology, University of Michigan.
