Dandome "contentillo" y justificando mi pereza, he de contarles que estos 5 meses no fueron del todo vanos, y no abandone completamente mi fiel y mas que justa causa, al difundir y vender R como un buen amante...para los solitarios habidos de software libre y sencillo.
Así pues, durante el " II Seminario de Biodiversidad y Conservación de Especies amenazadas, humedales, ecosistemas críticos", organizado por el "Grupo de Investigación y Estudios en Ciencias Biológicas" (Con el cual estoy muy agradecido, especialmente con Fernando Cediel, quien tuvo la "loca e irracional" idea de invitarme y con Natalia Rey, quien estuvo por ahí pendiente de los materiales necesarios para el curso (incluido el refrigerio! :P )) , tuve la oportunidad de dictar mi primer curso de manejo y estadística básica en R.
El curso tuvo una duración de 3 días (5, 6 y 7 de Septiembre, 2 horas cada día) y una asistencia de aproximadamente 28 estudiantes de distintos lugares del país; En un principio estuve un poquitin asustado con la idea, pero por la gracia de la buena onda, todo salio muy bien!
Gran parte del material que utilice en el curso provino de "R para Chibchombianos", y otra parte es complementaria a cabos que tenia sueltos (por ejemplo, la terminación de los post de exploración de los datos, previa a un análisis estadístico, basado en el articulo de Alain Zuur y colaboradores), para de esta forma reunir toda la información mostrada en el blog y otra adicional, compilada en un solo lugar.
De este modo, en este post, quiero compartirles via Google Docs (POR CIERTO, VIVA GOOGLE DOCS!!) todas las diapositivas que utilice en el curso, con comandos, información y todo el resto de cosas adicionales.
Las diapositivas contienen aveces imágenes escogidas por mi, de pintores que me gustan y fotógrafos que me gusta poner para hacer las presentaciones mas "firmes" y para compartir Arte, Ciencia, Software y vida juntos!!
La presentación esta organizada de la siguiente manera:
*Primer día: INTRODUCCIÓN AL PROGRAMA ESTADÍSTICO R
-Presentación
-Instalación de R
-Paquetes
-Repositorios
-Argumentos y comandos
-Objetos
-Vectores
-Matrices
-Arrays
-Data Frames
-Listas
-Factores
-Entrada y manejo de datos (incluido desde excel)
-Exportando datos y gráficos desde R
*Segundo día: ESTADÍSTICA DESCRIPTIVA Y UNIVARIADA EN R
-Frecuencias e Histogramas
-Tendencia central
-Normalidad
-Homogeneidad de varianzas
-Outliers
-Mas exploración de datos (recomendadisimo!!)
-T y T pareada
-U y U pareada
-Chi cuadrado
-Test de proporciones
-Anova de un factor
*Tercer día: ESTADÍSTICA MULTIVARIADA EN R
-Anova de dos factores
-Manova
-Kruskal-Wallis & “adonis”
-Clusters (bastante de esto!)
-Análisis discriminante
-Análisis de correspondencia canónica
-Componentes principales (PCA)
-Escalamiento multidimensional simple
-Perspectivas y ejercicios de R
ESPERO LES SEA DE GRAN AYUDA ESTAS DIAPOSITIVAS, EN ESPECIAL PARA MUCHAS COSAS QUE AÚN NO HABÍA COMPARTIDO CON UDS.
TODAVÍA NO SE SI ALCANCE A HACER EL ULTIMO POST DEL AÑO O SI ME TOCA QUEDÁRSELOS DEBIENDO PARA EL OTRO AÑO!
DE TODOS MODOS Y POR SI LOS TABANOS (LAS MOSCAS), LES AGRADEZCO MUCHO A LOS SEGUIDORES DEL POST, A LAS PERSONAS QUE ME LEEN POR AHÍ Y ESPERO QUE EL OTRO AÑO SEA REALMENTE PRODUCTIVO PARA TODOS EN CONJUNTO Y EN EN INDIVIDUAL!!
R para Chibchombianos-como utilizar R de la forma más desmenuzada posible
martes, 4 de diciembre de 2012
viernes, 23 de noviembre de 2012
El futuro próximo y/o lejano de R para Chibchombianos, después de dos años de cambiar pañales!
Después de 5 meses de inactividad en el blog, vuelve "R para Chibchombianos" la semana próxima con la intención de finalizar un ciclo que desde hace rato estaba sin concluir (me refiero al manejo básico de R, fundamentos de estadística descriptiva, univariada y multivariada).
Esto simplemente con la intención de dejar atrás el manejo básico de R, para adentrarnos en temas un poco mas avanzados, ya que el próximo 16 de Diciembre el blog cumplirá dos años de haberse creado y creo que el primer objetivo que era el entendimiento de las cosas básicas de R esta completo (o por lo menos lo estará la próxima semana).
Seguramente muchas cosas habrán hecho falta, pero siento que es imposible llegar a cubrirlas todas, solo espero que haya sido de gran ayuda el blog para las personas que no tenían ni idea de las ventajas de utilizar R y que haya sido una muy buena y muy fácil forma de empezar a entender R.
Que vendrá para la posteridad?
Bien, habiendo completado esta parte, creo que ahora me dedicare a lo mio y empezare a abordar cosas mas avanzadas, tanto para Biología en particular como para el publico en general:
-Ecología
-Sistemática Filogenética y taxonomía
-Manejo de datos geográficos
-Morfometría geométrica !!!!
-modelos estadísticos en general y ecológicos
-Estadística no-frecuentista
-funciones, código y lenguaje de programación !!!!
La transición hacia estos nuevos objetivos del blog sera muy divertida y en ocasiones también sera un poquito tediosa, sobre todo lo relacionado con programación (puesto que es de lo que menos se), pero sera quizás una muy buena forma de aprender muchas cosas juntos.
Espero que en los años venideros sea muy productivo el blog (mas de lo que pudo ser estos dos primeros años) y que les guste todas las cosas que van a ver.
NO OLVIDEN LA PRÓXIMA SEMANA "OJEAR" EL ULTIMO POST DE ESTE PRIMER CICLO, SERÁ INTERESANTE (BUENO, ESO LO DIRÁN USTEDES, PARA MI LO ES!), MUY CONCISO Y CONTENDRÁ MUCHA INFORMACIÓN (YA DISPONIBLE EN EL BLOG Y NUEVA INFORMACIÓN) EN UN ESPACIO REDUCIDO, A MANERA DE RESUMEN DE TODA LA PRIMER PARTE DE "R PARA CHIBCHOMBIANOS"!!
"firme"-mente
BROSSY
sábado, 16 de junio de 2012
1/1000 Maneras de Graficar un Mapa en R
Que tal chibchombianos!
Así como existen 1000 maneras de morir tambien podemos hablar de que existen 1000 formas de hacer todo en R y en este caso graficar un mapa no es la excepción. Aunque personalmente me tomo tiempo encontrar una forma rápida y sencilla para lograrlo, aqui les ahorro esfuerzo y sufrimiento!...
Simplemente deben seguir las siguientes pasos para ser felices:
1. Descargar el shapefile.
El shapefile es un formato para sistemas de información geográfica. De acuerdo con la Wiki es un formato vectorial de almacenamiento digital donde se
guarda la localización de los elementos geográficos y los atributos
asociados a ellos. Si bien la definición es básica, es increible la cantidad de información descriptiva de un mapa que se puede almacenar en un shapefile.
Existe un shapefile disponible casi que para cualquier región y es posible encontrar varios sitios en la web que permiten acceder gratuitamente a sus shapefiles, en este caso les recomiendo la base de datos del GADM a la cual podrán acceder con el siguiente link http://www.gadm.org/
Otros links de interés
Para el ejemplo utilizaremos la base de datos del GADM, asi que vayan a la sección de descarga y seleccionen el mapa de Uzbekistan (de paso se enteran que existe este país) y descarguen el shapefile respectivo.
1a. Instalen y carguen en R las bibliotecas maptools, sp y RColorBrewer
Para instalar los paquetes:
> install.packages("maptools",lib="/directorio")
Para instalar los paquetes:
> install.packages("maptools",lib="/directorio")
Para cargar:
> library(maptools)
# Provee herramientas para leer y manejar objetos espaciales
# Provee herramientas para leer y manejar objetos espaciales
> library(sp)
#Provee métodos para manejar datos espaciales
#Provee métodos para manejar datos espaciales
> library(RcolorBrewer)
# Permite acceder a paletas de colores para cartografía
# Permite acceder a paletas de colores para cartografía
1b. Ahora lean el shapefile y conviertanlo en un vector, en este caso lo llamé ohsi
> ohsi<-readShapeSpatial("UZB_adm1.shp")
#Corresponde al mapa politico de Uzbekistan
#Corresponde al mapa politico de Uzbekistan
1c. Comprueben que se cargo el mapa
> plot(ohsi)
1d. Pueden incluso revisar la estructura interna del vector
> str(ohsi)
> ohsi@data
1.e. Asignen el vector como una variable
> mapa=ohsi
2. El siguiente paso es construir una tabla que contenga en una columna las localidades "LOCAL" y en la siguiente columna algún criterio para la selección de las localidades, en este caso el criterio fue "lugares que quiero visitar en Uzbekistan" y la llame "VISITA" :P. Es importante guardar la tabla en formato csv o delimitado por espacios.
Aqui tienen la tabla
"lOCAL" "VISITA"
"Andijon" "NA"
"Bukhoro" 1
"Ferghana" "NA"
"Jizzakh" "NA"
"Karakalpakstan" "NA"
"Kashkadarya" "NA"
"Khorezm" "NA"
"Namangan" "NA"
"Navoi" "NA"
"Samarkand" 1
"Sirdaryo" "NA"
"Surkhandarya" "NA"
"Tashkent City" "NA"
"Tashkent" 1
"Andijon" "NA"
"Bukhoro" 1
"Ferghana" "NA"
"Jizzakh" "NA"
"Karakalpakstan" "NA"
"Kashkadarya" "NA"
"Khorezm" "NA"
"Namangan" "NA"
"Navoi" "NA"
"Samarkand" 1
"Sirdaryo" "NA"
"Surkhandarya" "NA"
"Tashkent City" "NA"
"Tashkent" 1
La codificación corresponde a: NA=No aplica y 1=lugares que quiero visitar
2a. Ahora deben abrir el archivo de la tabla y crear otro vector el cual llamé "yque"
> yque<-read.csv("tabla.csv",header =TRUE, sep = "\t")
2b. Revisen los datos de la columna visita en la tabla. Deberian ser visibles en la consola.
> yque$VISITA
2c. Reasignen la información de vector "yque" variable "VISITA" a "mapa"
> mapa@data=data.frame(yque$VISITA)
> mapa
3. Aunque hay múltiples paquetes para asignar colores en un mapa, en este caso me base en la paleta de colores de colorBrewer. El color pueden seleccionarlo guiándose con el siguiente link http://colorbrewer2.org/
Ahora grafiquen el mapa con las localidades seleccionadas en "VISITA" y el color lo deben especificar en col.regions=brewer.pal() como sigue a continuación:
> spplot(mapa,c("yque.VISITA"),col.regions=brewer.pal(3, "RdYlBu"), scales=list(draw = TRUE))
Listo!! ahora tenemos el mapa de Uzbekistan y las áreas seleccionadas corresponden a las localidades que me gustarian visitar :P.
otro ejemplo
Fácil! Ahora lo pueden modicar y adaptar a sus datos. Esta es solo una forma sencilla, les recuerdo hay cientos de opciones por explorar!
Buena energia!
jueves, 22 de marzo de 2012
Errores comunes en estadística: Exploración de datos (I) en R: Outliers, Homogeneidad y Normalidad
En marzo de 2010, Alain Zuur, Elena Leno y Chris Elphick, publicaron en la revista gratuita "Methods in ecology and evolution", un artículo supremamente interesante, muy sencillo, pero de gran importancia.
En este presentan un protocolo sencillo que consta de 7 pasos en la exploración de datos para evitar errores comunes en estadística que son pasados por alto cuando se analiza un set de datos.
En muchas ocasiones este tipo de desaciertos pueden llevarnos a cometer errores tipo I o tipo II, y es "altamente probable" que esto afecte nuestros resultados y por consiguiente nuestras conclusiones.
Como Chibchombia es un país lleno de "falsos positivos", debemos suponer que no importa si aparecen uno o dos falsos positivos mas diarios; pero no es así!!. Nos importan de sobremanera.
Entonces no nos conviene pues a nosotros cometer ni un falso positivo mas, así sea en estadística.
En este primer post sobre exploración de datos, empezaré describiendo y mostrando el camino para realizar en R los tres primeros pasos del protocolo propuesto por Zuur y colaboradores, estos tres primero pasos son:
1. Búsqueda de Outliers en las variables del set de datos.
2. Prueba de homogeneidad de varianzas.
3. Prueba de distribución normal de los datos.
En 2 posts posteriores describiré los 4 pasos restantes del protocolo. Por ahora espero les guste este post y aun mas les sea muy útil.
Lo "firme*" de este protocolo de exploración de datos, además de su amplia utilidad, es que los resultados que se obtienen son gráficos, y se utilizan varios tipos de diagramas disponibles en paquetes de R.
(En la medida de lo posible también incluiré algunos caminos de exploración de datos de forma no gráfica disponibles en R; En el artículo solo se presentan los resultados de modo gráfico, pero aveces es útil obtener resultados numéricos como tal (p. ej. en pruebas de homogeneidad de varianza y normalidad)).
*FIRME= en el lenguaje "ñeristico", significa, chevere, bacano, bueno, cool, etc!
Sin mas "formalismos", empecemos:
1. Búsqueda de Outliers en las variables del set de datos.
El objeto primordial de este blog no es explicar las implicaciones estadísticas de las técnicas utilizadas, pero debo recordar que cada investigador debe conocer las consecuencias y el impacto que tiene la presencia de outliers al utilizar una técnica estadística en particular como NDMS o Modelos lineales de Poisson.
Estoy limitado por tanto a mencionar que un outlier es una observación/dato que tiene un valor relativo mucho mayor o menor comparado con el resto de observaciones del set de datos; Y por supuesto a explicar a cabalidad algunos de los caminos para realizar la búsqueda de outliers en R, utilizados por los autores en el artículo.
En el artículo se utiliza una base de datos grande de las medidas presentes en Ammodramus caudacutus, esta matriz de datos se encuentra disponible en este link.
Pero para entender mejor el procedimiento (y para utilizar una base de datos más pequeña de sólo 190 datos, puesto que la base de datos del artículo es de más de 1000 observaciones) construiremos nosotros nuestra propia base de datos tomando la medida de longitud estandar (LTS) de Chaetostoma leucomelas.
Para esto sólo tenemos que escribir un archivo de texto que contenga lo siguiente, con el nombre de "chae.txt": (se debe crear el archivo en una sola columna por lo tanto, los siguientes datos todos pertenecen a una sola variable "lts"; Escribirlos en el .txt como una sola columna hacia abajo!!!)
De este modo, teniendo nuestra base de datos en un archivo .txt, procedemos a llamarla en R (para entender más sobre entrada y manejo de datos en R ir a Entrada y manejo de datos en R (guía básica).
> c_leucomelas <- read.table(file = "chae.txt", header = TRUE)
# Llamamos al objeto creado en R "c_leucomelas", utilizamos el comando read.table puesto que leeremos un archivo .txt y le decimos que la primera fila es el nombre de la variable "lts", utilizando header=TRUE.
Para la búsqueda de outliers gráficamente utilizaremos los dos métodos usados por los autores en el artículo:
1. Gráfica de caja y bigotes (boxplot)
2. Gráfica de cleveland (cleveland dotplot)
> boxplot(c_leucomelas$lts, ylab = "longitud estandar (mm)")
# Realizamos nuestro boxplot, llamando la variable "lts" del objeto "c_leucomelas" con el argumento "c_leucomelas$lts" y a nuestro eje Y le ponemos la etiqueta "longitud estandar (mm)"
> dotchart(c_leucomelas$lts, xlab = "longitud estandar (mm)", ylab = "Orden de los datos")
#Realizamos nuestro cleveland dotplot, pero en este caso le ponemos como etiqueta al eje x "longitud estandar (mm)" y al eje Y "orden de los datos".
Básicamente la diferencia entre los dos tipos de gráfico, es que el cleveland dotplot, nos especifica el orden de los datos, lo cual nos ayudaria a detectar particularmente cuales datos son outliers.
Ahora para ir un poco mas alla y obtener gráficos bien bonitos y funcionales, utilizaremos otra línea de comandos que nos grafique el cleveland dotplot y el boxplot en uno sólo, y utilizaremos el paquete "lattice" junto con la base de datos del artículo, para realizar un cleveland dotplot multiple con todas las variables!
> par(mfrow= c (1,2), mar = c(5,4,2,1))
# Se utiliza el argumento "par" para indicar que se incluiran varios gráficos en uno solo; el argumento "mfrow" indica el número de gráficos y su disposición espacial en el gráfico total, en nuestro ejemplo ponemos "= c(1,2)" lo que indica que habrá 1 fila y dos columnas, si quisieramos poner el gráfico de manera vertical tendriamos que darle "= c(2,1)" lo cual indica que los dos gráficos que vamos a poner van uno encima del otro, es decir 2 filas y una columna.
La segunda parte del comando indica el tamaño de los margenes de la figura total, en el primer margen que es el de abajo le damos un valor de 5, en el segundo que es el de la izquierda le damos un valor de 4, en el tercero que es el de arriba le damos un valor de 2 y finalmente en el cuarto que corresponde al margen de la derecha le damos un valor de 1.
> boxplot(c_leucomelas$lts, ylab = "lts (mm)")
# Realizo el gráfico de cajas y bigotes
> dotchart(c_leucomelas$lts, xlab = "lts (mm)",ylab = "Orden de los datos")
#Realizo el gráfico de Cleveland
>install.packages ("lattice")
# Instalo el paquete "lattice" en R. Hay que tener cuidado acerca de la version de R que estemos utilizando, puesto que lattice esta actualizandose muy seguido y si tenemos una version de R anterior, es necesario descargar una version anterior del paquete lattice desde la pagina y cargarla a nuestra versión de R.
Para ver todas las versiones de lattice pueden ir al siguiente link.
>library(lattice)
#Cargo el paquete "lattice" en mi área de trabajo de R.
>Sparrows<- read.table(file = "SparrowsElphick.txt", header = TRUE)
# Llamo la matriz de datos a R, con el nombre "Sparrows".
>super<- cbind(Sparrows$wingcrd, Sparrows$tarsus, Sparrows$head,Sparrows$culmen, Sparrows$nalospi,Sparrows$wt)
# Creo un compilado de las variables que utilizaremos para hacer nuestro cleveland dotplot y que extraeremos de la base de datos de los autores, a este compilado de variables le llamaremos "super".
>colnames(super) <- c("long. ala", "long. tarso", "long. cabeza", "long. mulmen", "long. total", "peso")
# Creamos un vector con los nombres (colnames) de las variables ordenadas de nuestro compilado (super) creado anteriormente.
>dotplot(as.matrix(super), groups = FALSE,
strip = strip.custom(bg = 'green',
par.strip.text = list(cex = 0.8)),
scales = list(x = list(relation = "free"),
y = list(relation = "free"),
draw = FALSE),
col = 1, cex = 0.5, pch = 16,
xlab = "valor de la variable",
ylab = "orden de los datos en el archivo texto")
# Para empezar, le indicamos a R el tipo de gráfico que vamos a realizar (dotplot), y le decimos que organice todos los gráficos en una matriz de nuestras variables (super). Con el argumento "groups= FALSE" le indicamos a R que queremos todas las variables en gráficos diferentes, es decir cada una con un eje X y un eje Y independientes, sí le diéramos "groups=TRUE" nos graficaria todas las variables con un mismo eje X y un mismo eje Y, osea que nos las agruparía. Después lo que vamos a hacer es modificar la apariencia del gráfico con los argumentos de "strip", el primero es el color del fondo de los recuadros de los nombre de las variables, (bg=background= fondo), yo utilice el color verde "bg = 'green'" pero uds podrian utilizar cualquier color, como el amarillo, rojo o azul, recuerden que se debe escribir en ingles; Después se definen las escalas de la lista de valores en los ejes X y Y de cada variable, en nuestro gráfico nosotros dejaremos las variables con una relación libre "free", puesto que no las vamos a utilizar, si quisiéramos que fueran puestas en el diagrama, tendríamos que darle " draw=TRUE", el problema es que como tenemos una gran cantidad de datos los nombres de cada uno de los datos del eje Y se superpondrán en la gráfica y no se entenderá nada. Finalmente lo que hacemos es modificar los puntos que utilizaremos dentro de las gráficas, con el argumento "col=1" le decimos a R el numero de color determinado que queremos utilizar en nuestro gráfico para los puntos del gráfico, nosotros utilizamos en color negro que es codificado con el numero 1, si lo quisiéramos de color rojo tendríamos que darle el numero 2 y si quisiéramos de color verde le daríamos 3 y asi sucesivamente, despues lo que hacemos sera definir el tamaño de los puntos que utilicemos en la grafica, nosotros utilizaremos un tamaño de 0.5 "cex= 0.5" pero eso depende del gusto de cada quien, tambien podemos cambiar el simbolo que queremos utilizar, por ejemplo pueden ponerse asteriscos, estrellas, triangulos o simplemente puntos como lo hicimos aqui "pch=16", si hubiésemos querido poner estrellas en vez de puntos debíamos darle "pch=11", se pueden encontrar todos los valores para cada simbolo que utilicemos en la pagina de quick R. Ya lo ultimo es poner las etiquetas de los dos ejes X y Y (xlab = "valor de la variable", ylab = "orden de los datos en el archivo texto") .
Teniendo en cuenta que son dos funciones en una sola, lo primero que hacemos es darle un nombre a la función 1, por lo tanto llamaremos nuestra función 1 "buscar_outliers", esta función sera aplicada para la matriz "X" (buscar_outliers<- function (x)) que puede ser cualquier matriz, y lo que hace es aplicar la función 2 a todas las columnas de la matriz (x), para esto se utiliza el argumento "return (apply (x,2,buscar))" que indica lo que queremos obtener como resultado (return) de la función que en este caso son las posiciones de los outliers en cada una de las variables, por ello le indicamos que busque los outliers de las columnas/variables con el numero "2", si quisieramos buscar los outliers del las filas tendriamos que darle el numero "1", finalmente la palabra buscar es el nombre de la función 2 e indica que se aplique la funcion numero 2 a cada una de las columnas de la matriz x.
# De este modo, lo único que tenemos que dar es el nombre de nuestra función creada y especificar a cual matriz queremos aplicarla que en nuestro caso es "Sparrows", y nuestros resultados nos dirán cual es el dato en cada una de las variables que posee un valor outlier (por ejemplo para "wingcrd" el dato numero 536 debe ser considerado un outlier):
>buscar_outliers (Sparrows)
Esto que acabo de presentar es un pequeñisimo ejemplo de como construir funciones, pero en la vida real no es necesario construir funciones para encontrar outliers (aunque es necesario construir muchísimas funciones para nuestras necesidades particulares), es importante empezar con algo bastante sencillo.
De hecho, existen bastantes paquetes disponibles en R para buscar outliers, PARA MI uno muy bueno y fácil de utilizar es el paquete "outliers", daré un par de ejemplos de las funciones disponibles es este paquete aplicándolas a nuestra matriz.
>install.packages ("outliers")
#Instalamos el paquete "outliers"
> library (outliers)
#Cargamos el paquete en nuestra área de trabajo.
>chisq.out.test(Sparrows$wingcrd)
# Realizo una búsqueda de outlier por medio del método de chi-cuadrado, este nos indicara el valor de la variable que debe ser tomada como un outlier en nuestra variable "wingcrd" de nuestra matriz "Sparrows"
(Para mas info acerca de este método...)
El método gráfico que utilizaron los autores fue el método de "cajas y bigotes condicionales" (conditional boxplots) y este es el que voy a mostrar en el post ademas de mostrar una pequeñísima linea de comandos para realizar el test de Levene para homogeneidad de varianzas.
Para esta prueba utilizaremos la base de datos de los autores, esta contiene las tasas de consumo de alimento de una ave migradora (Limosa haemastica) y se quiere saber si el consumo de alimento varia entre sexos, periodos de tiempo o por una combinación de estas dos variables. Por lo tanto se debe asumir que la (i) variacion en las obseraciones delos sexos es similar, (ii) que la variación de las observaciones en diferentes periodos de tiempo es similar y (iii) que la variación entre periodos de tiempo y sexos es similar.
> library (lattice)
#Se carga el paquete "lattice" a nuestra área de trabajo.
>limosa <- read.table(file="Godwits.txt", header=TRUE)
#Llamamos a R la base de datos para hacer nuestro análisis, y le ponemos por nombre "limosa", el archivo en el que están guardado los datos fue llamado por los autores "Godwits.txt" y le damos "header=TRUE" para indicarle que la primera fila son los nombres de las variables.
>limosa$fSEX <- factor(limosa$SEX, levels = c(0, 1, 2),labels = c("No", "Hembras", "Machos"))
# Debido a que los datos están organizados de forma codificada, debemos crear un vector de nombres para los estados en la tabla; es decir, en el archivo para la variable "SEX" el estado 0 significa que tiene sexo desconocido, 1 que es una hembra y 2 que el individuo es un macho. Por ello se le dice a R que construya un factor "SEX" de la base de datos "limosa" (limosa$fSEX) a partir de los estados 0,1 y 2 de la variable "SEX" de la base "limosa" (limosa$SEX, levels = c(0, 1, 2)) y que etiquete cada uno de los estados respectivamente con los nombres "No", "Hembras" y "Machos" (labels = c("No", "Hembras", "Machos").
>limosa$fPERIOD <- factor(limosa$PERIOD, levels = c(0, 1, 2),labels = c("Verano", "Pre-migracion", "Invierno"))
#Se hace lo mismo que en el paso anterior, solo que esta vez lo que codificaremos es los tres periodos de tiempo con los nombres "Verano", "Pre-migracion" e "Invierno".
Finalmente lo que hacemos es la parte de edición del gráfico para que quede mas bonito o quede como nosotros queramos, para ello la linea de argumentos que utilizamos es la de ajuste de los parámetros del gráfico especificando con una lista los parámetros que queramos editar (par.settings =list(.....)). En nuestro gráfico solo editaremos 3 parámetros del gráfico (aunque uno puede editar muchísimos parámetros y solo haría falta adicionarlos a la lista), el primer parámetro es el color del rectángulo de las cajas, al cual le dimos un color verde que esta codificado con el numero "3", si quisiéramos aplicar el color rojo tendríamos que darle el numero "2" (box.rectangle = list(col = 3)) . El segundo parámetro que modificaremos es el color de el bigote de la caja (en ingles, la sombrilla de la caja), le damos el color numero "2" que es rojo (box.umbrella = list(col = 2)) . El tercer parámetro que vamos a editar es el del color de los puntos extremos (outliers) de las variables y el tamaño de estos en el gráfico, entonces le damos un tamaño de punto de 1.5 y elegimos el color negro para estos que esta codificado con el numero "1" (plot.symbol = list(cex = 1.5, col = 1)). Finalmente lo que hacemos es determinar la escala de los ejes pero en este caso el decimos que la relación de escalas entre los ejes es la misma (scales = list(x = list(relation = "same"), y = list(relation = "same"))) . Las etiquetas (labels) de los ejes fueron nombradas con el comando (xlab = "Periodo de migracion", ylab = "Tasa de consumo de alimento").
>grafico_homogeneidad
# Visualizamos nuestro gráfico.
Ahora bien, al igual que para encontrar outliers existen métodos no gráficos, para revisar sí tenemos homogeneidad de varianza en nuestros datos también hay muchos test disponibles en R, como la prueba de F que es la que utilizare para explicar su desarrollo en R, además del test de Levene entre otros.
> hembras<-subset(limosa, SEX==1, select = c(mgconsumed),header =T)
#Primero hacemos un substet que contenga la cantidad de alimento consumido por las hembras, utilizando la base de datos "limosa", el estado para la variable "SEX" utilizado para las hembras es el numero 1 (SEX==1), y seleccionamos los valores de cantidad de alimento consumido (select = c(mgconsumed)), correspondiente a los individuos que son hembras. Por ultimo le decimos que la primera fila de "limosa" es el nombre de las variables y también del nuevo subset (header =T).
> machos<-subset(limosa, SEX==2, select = c(mgconsumed),header =T)
#Hacemos lo propio, creando también el subset de los machos para comparar las varianzas.
>var.test(hembras$mgconsumed,machos$mgconsumed)
# Realizo mi test de homogeneidad de varianzas para los dos substes antes creados y obtenog mis resultados:
3. Prueba de distribución normal de los datos.
Finalmente, para terminar este post tan largo, pero el cual creo que es muy útil (además es demasiado largo porque me esforcé muchísimo por explicar cada detalle minuciosamente, así que espero todo haya quedado muy claro), explicare la vía gráfica (histograma de frecuencias) para hacer la prueba de distribución normal en R, además de explicar una vía no gráfica entre una infinidad de formas de hacerlo en R.
>pericos <- read.table(file="SparrowsElphick.txt", header=TRUE)
#Para empezar lo primero que vamos a hacer es llamar nuestros datos (los datos de los autores, y crearemos un objeto en R llamado "pericos" (no son pericos, pero que hijuemadre, pongamole pericos), y le decimos que la primer fila de la tabla son los nombres de las variables con (header=T).
>pericos$fMonth<-factor(pericos$Month,
levels = c(5, 6, 7, 8, 9, 10),
labels = c("Mayo", "Junio", "Julio", "Agosto",
"Sept.", "Oct."))
#Despues creo un subset de datos, y para ello lo primero que debo hacer es escoger un factor que me separe algunas observaciones de otras, por lo tanto escojo la variable "mes" (month) del 5,6,7,8,9 y 10 que es como están codificados en la matriz y corresponden a los meses de mayo hasta Octubre) set de datos original, y extraigo las observaciones correspondientes de estos 6 meses.
>pericos$I1 <- pericos$fMonth =="Junio" |
pericos$fMonth =="Julio" |
pericos$fMonth =="Agosto" |
pericos$fMonth =="Mayo" |
pericos$fMonth =="Oct." |
pericos$fMonth =="Sept."
# Ahora creo un objeto en R con las observaciones obtenidas anteriormente, pero esta vez, cuando los mese ya han sido etiquetados y no están codificados numéricamente.
>hist(pericos$wt[pericos$I1],
xlab = "Peso (g)", breaks = 100,
main = "", ylab = "Frecuencia")
# Finalmente realizamos el histograma de frecuencias para visualizar si existe o no normalidad en los datos, utilizando los datos de peso de la matriz original (pericos$wt) (wt=weight=peso) de las observaciones correspondientes al factor antes creado que corresponde a algunos meses ([pericos$I1]), le ponemos nombre a los ejes con (xlab = "Peso (g) y "ylab = "Frecuencia"), y le decimos que el histogramas de frecuencias haga 100 rangos (breaks = 100), es decir 100 barritas de un peso en especifico cada una, pero podriamos decirle que haga el numero que nosotros deseemos.
De vez en cuando se deben dividir los datos, para así definir que parte esta cargando la mayor variabilidad o como se comportan los datos, por ejemplo en meses diferentes, para ello podemos hacer un histograma de frecuencias, para así revisar la normalidad de los datos de un mes en especifico.
>library(lattice)
#Cargo el paquete "lattice", para hacer graficos multiples.
>histogram( ~ wt | fMonth, type = "count",
xlab = "Peso (g)",
ylab = "Frecuencia",
nint=100,layout=c(3,2),
strip.left = strip.custom(bg = 'yellow'),
strip = F,
col.line = "red", col = "red",
scales = list(x = list(relation = "same"),
y = list(relation = "same"),
draw = TRUE),
subset = fMonth =="Mayo" |fMonth =="Junio" | fMonth == "Julio" |fMonth == "Agosto" |fMonth == "Sept." |fMonth == "Oct.",
data = pericos)
#Realizo mi histograma de frecuencias con los datos de peso (wt) de mi vector fMonth, y el tipo de histograma a realizar es del tipo "count", le pongo el nombre de los ejes de cada variable, con el argumento (nint=100), le digo que divida cada histograma en 100 rangos, con (layout=c(3,2), le digo que distribuya cada uno de los histogramas en 3 columnas y dos filas, los siguientes argumentos son para cambiar los colores de los gráficos las escalas de cada uno de los ejes para dejar las mismas para todos los gráficos; por ultimo especificamos el subset de datos que vamos a utilizar.
Y para realizar un aprueba de normalidad rápidamente en R, hay muchos caminos pero el que generalmente he usado yo es el test de Shapiro, disponible en algunos de mis posts anteriores.
uffffff
finalmente eso era todo lo de esta primera parte de exploración de datos para evitar errores estadísticos, y para los próximos dos posts les tendré los pasos restantes de todo el protocolo.
Proximo post:
Errores comunes en estadística: Exploración de datos (II) en R: Cantidad de ceros en los datos y Colinealidad entre las variables.
En este presentan un protocolo sencillo que consta de 7 pasos en la exploración de datos para evitar errores comunes en estadística que son pasados por alto cuando se analiza un set de datos.
En muchas ocasiones este tipo de desaciertos pueden llevarnos a cometer errores tipo I o tipo II, y es "altamente probable" que esto afecte nuestros resultados y por consiguiente nuestras conclusiones.
Como Chibchombia es un país lleno de "falsos positivos", debemos suponer que no importa si aparecen uno o dos falsos positivos mas diarios; pero no es así!!. Nos importan de sobremanera.
Entonces no nos conviene pues a nosotros cometer ni un falso positivo mas, así sea en estadística.
En este primer post sobre exploración de datos, empezaré describiendo y mostrando el camino para realizar en R los tres primeros pasos del protocolo propuesto por Zuur y colaboradores, estos tres primero pasos son:
1. Búsqueda de Outliers en las variables del set de datos.
2. Prueba de homogeneidad de varianzas.
3. Prueba de distribución normal de los datos.
En 2 posts posteriores describiré los 4 pasos restantes del protocolo. Por ahora espero les guste este post y aun mas les sea muy útil.
Lo "firme*" de este protocolo de exploración de datos, además de su amplia utilidad, es que los resultados que se obtienen son gráficos, y se utilizan varios tipos de diagramas disponibles en paquetes de R.
(En la medida de lo posible también incluiré algunos caminos de exploración de datos de forma no gráfica disponibles en R; En el artículo solo se presentan los resultados de modo gráfico, pero aveces es útil obtener resultados numéricos como tal (p. ej. en pruebas de homogeneidad de varianza y normalidad)).
*FIRME= en el lenguaje "ñeristico", significa, chevere, bacano, bueno, cool, etc!
Sin mas "formalismos", empecemos:
1. Búsqueda de Outliers en las variables del set de datos.
El objeto primordial de este blog no es explicar las implicaciones estadísticas de las técnicas utilizadas, pero debo recordar que cada investigador debe conocer las consecuencias y el impacto que tiene la presencia de outliers al utilizar una técnica estadística en particular como NDMS o Modelos lineales de Poisson.
Estoy limitado por tanto a mencionar que un outlier es una observación/dato que tiene un valor relativo mucho mayor o menor comparado con el resto de observaciones del set de datos; Y por supuesto a explicar a cabalidad algunos de los caminos para realizar la búsqueda de outliers en R, utilizados por los autores en el artículo.
En el artículo se utiliza una base de datos grande de las medidas presentes en Ammodramus caudacutus, esta matriz de datos se encuentra disponible en este link.
Pero para entender mejor el procedimiento (y para utilizar una base de datos más pequeña de sólo 190 datos, puesto que la base de datos del artículo es de más de 1000 observaciones) construiremos nosotros nuestra propia base de datos tomando la medida de longitud estandar (LTS) de Chaetostoma leucomelas.
Para esto sólo tenemos que escribir un archivo de texto que contenga lo siguiente, con el nombre de "chae.txt": (se debe crear el archivo en una sola columna por lo tanto, los siguientes datos todos pertenecen a una sola variable "lts"; Escribirlos en el .txt como una sola columna hacia abajo!!!)
| lts | 55.08 | 75.52 | 82.78 | 63.66 | 46.26 | 32.26 | 24.88 |
| 74.20 | 31.76 | 76.88 | 72.00 | 47.20 | 43.80 | 31.58 | 21.76 |
| 69.94 | 26.24 | 68.10 | 89.92 | 62.46 | 48.14 | 27.62 | 23.06 |
| 82.48 | 29.98 | 87.60 | 73.20 | 99.28 | 57.00 | 31.72 | 22.02 |
| 85.44 | 28.08 | 70.60 | 72.34 | 111.96 | 50.26 | 33.48 | 22.94 |
| 90.76 | 33.42 | 65.84 | 69.78 | 90.24 | 58.22 | 29.10 | 21.50 |
| 90.32 | 38.36 | 78.04 | 99.72 | 74.58 | 49.90 | 29.56 | 23.28 |
| 60.12 | 23.02 | 70.62 | 83.28 | 90.02 | 43.26 | 28.06 | 26.04 |
| 102.58 | 34.72 | 62.16 | 88.58 | 80.42 | 47.56 | 28.68 | 22.76 |
| 97.20 | 27.56 | 68.90 | 71.24 | 80.14 | 35.06 | 26.54 | 23.16 |
| 98.42 | 29.82 | 82.66 | 87.98 | 73.82 | 28.98 | 26.04 | 23.76 |
| 88.86 | 30.92 | 59.02 | 86.46 | 71.88 | 30.92 | 25.58 | 25.96 |
| 76.58 | 32.54 | 73.54 | 81.94 | 71.70 | 31.54 | 29.24 | 22.18 |
| 85.15 | 26.52 | 53.14 | 74.18 | 70.22 | 33.62 | 25.02 | 23.80 |
| 75.22 | 31.00 | 44.78 | 74.19 | 64.18 | 32.50 | 26.56 | 20.24 |
| 82.72 | 25.94 | 48.24 | 71.38 | 66.00 | 34.32 | 25.38 | 21.12 |
| 69.74 | 21.72 | 38.56 | 65.99 | 65.16 | 32.68 | 25.32 | 19.10 |
| 72.64 | 25.12 | 58.94 | 69.09 | 61.56 | 31.34 | 26.08 | 18.14 |
| 67.06 | 22.96 | 47.32 | 70.58 | 62.34 | 29.34 | 25.56 | 60.22 |
| 66.89 | 16.92 | 29.68 | 65.82 | 54.88 | 31.08 | 23.54 | 95.10 |
| 67.52 | 20.48 | 72.24 | 71.18 | 54.90 | 31.56 | 23.38 | 71.70 |
| 66.88 | 22.82 | 99.00 | 63.04 | 53.76 | 31.16 | 28.62 | |
| 59.51 | 65.72 | 98.54 | 56.69 | 48.14 | 30.44 | 23.20 | |
| 58.52 | 92.40 | 95.74 | 57.12 | 57.86 | 31.48 | 26.82 |
De este modo, teniendo nuestra base de datos en un archivo .txt, procedemos a llamarla en R (para entender más sobre entrada y manejo de datos en R ir a Entrada y manejo de datos en R (guía básica).
> c_leucomelas <- read.table(file = "chae.txt", header = TRUE)
# Llamamos al objeto creado en R "c_leucomelas", utilizamos el comando read.table puesto que leeremos un archivo .txt y le decimos que la primera fila es el nombre de la variable "lts", utilizando header=TRUE.
Para la búsqueda de outliers gráficamente utilizaremos los dos métodos usados por los autores en el artículo:
1. Gráfica de caja y bigotes (boxplot)
2. Gráfica de cleveland (cleveland dotplot)
> boxplot(c_leucomelas$lts, ylab = "longitud estandar (mm)")
# Realizamos nuestro boxplot, llamando la variable "lts" del objeto "c_leucomelas" con el argumento "c_leucomelas$lts" y a nuestro eje Y le ponemos la etiqueta "longitud estandar (mm)"
> dotchart(c_leucomelas$lts, xlab = "longitud estandar (mm)", ylab = "Orden de los datos")
#Realizamos nuestro cleveland dotplot, pero en este caso le ponemos como etiqueta al eje x "longitud estandar (mm)" y al eje Y "orden de los datos".
Básicamente la diferencia entre los dos tipos de gráfico, es que el cleveland dotplot, nos especifica el orden de los datos, lo cual nos ayudaria a detectar particularmente cuales datos son outliers.
Ahora para ir un poco mas alla y obtener gráficos bien bonitos y funcionales, utilizaremos otra línea de comandos que nos grafique el cleveland dotplot y el boxplot en uno sólo, y utilizaremos el paquete "lattice" junto con la base de datos del artículo, para realizar un cleveland dotplot multiple con todas las variables!
> par(mfrow= c (1,2), mar = c(5,4,2,1))
# Se utiliza el argumento "par" para indicar que se incluiran varios gráficos en uno solo; el argumento "mfrow" indica el número de gráficos y su disposición espacial en el gráfico total, en nuestro ejemplo ponemos "= c(1,2)" lo que indica que habrá 1 fila y dos columnas, si quisieramos poner el gráfico de manera vertical tendriamos que darle "= c(2,1)" lo cual indica que los dos gráficos que vamos a poner van uno encima del otro, es decir 2 filas y una columna.
La segunda parte del comando indica el tamaño de los margenes de la figura total, en el primer margen que es el de abajo le damos un valor de 5, en el segundo que es el de la izquierda le damos un valor de 4, en el tercero que es el de arriba le damos un valor de 2 y finalmente en el cuarto que corresponde al margen de la derecha le damos un valor de 1.
> boxplot(c_leucomelas$lts, ylab = "lts (mm)")
# Realizo el gráfico de cajas y bigotes
> dotchart(c_leucomelas$lts, xlab = "lts (mm)",ylab = "Orden de los datos")
#Realizo el gráfico de Cleveland
>install.packages ("lattice")
# Instalo el paquete "lattice" en R. Hay que tener cuidado acerca de la version de R que estemos utilizando, puesto que lattice esta actualizandose muy seguido y si tenemos una version de R anterior, es necesario descargar una version anterior del paquete lattice desde la pagina y cargarla a nuestra versión de R.
Para ver todas las versiones de lattice pueden ir al siguiente link.
>library(lattice)
#Cargo el paquete "lattice" en mi área de trabajo de R.
>Sparrows<- read.table(file = "SparrowsElphick.txt", header = TRUE)
# Llamo la matriz de datos a R, con el nombre "Sparrows".
>super<- cbind(Sparrows$wingcrd, Sparrows$tarsus, Sparrows$head,Sparrows$culmen, Sparrows$nalospi,Sparrows$wt)
# Creo un compilado de las variables que utilizaremos para hacer nuestro cleveland dotplot y que extraeremos de la base de datos de los autores, a este compilado de variables le llamaremos "super".
>colnames(super) <- c("long. ala", "long. tarso", "long. cabeza", "long. mulmen", "long. total", "peso")
# Creamos un vector con los nombres (colnames) de las variables ordenadas de nuestro compilado (super) creado anteriormente.
>dotplot(as.matrix(super), groups = FALSE,
strip = strip.custom(bg = 'green',
par.strip.text = list(cex = 0.8)),
scales = list(x = list(relation = "free"),
y = list(relation = "free"),
draw = FALSE),
col = 1, cex = 0.5, pch = 16,
xlab = "valor de la variable",
ylab = "orden de los datos en el archivo texto")
# Para empezar, le indicamos a R el tipo de gráfico que vamos a realizar (dotplot), y le decimos que organice todos los gráficos en una matriz de nuestras variables (super). Con el argumento "groups= FALSE" le indicamos a R que queremos todas las variables en gráficos diferentes, es decir cada una con un eje X y un eje Y independientes, sí le diéramos "groups=TRUE" nos graficaria todas las variables con un mismo eje X y un mismo eje Y, osea que nos las agruparía. Después lo que vamos a hacer es modificar la apariencia del gráfico con los argumentos de "strip", el primero es el color del fondo de los recuadros de los nombre de las variables, (bg=background= fondo), yo utilice el color verde "bg = 'green'" pero uds podrian utilizar cualquier color, como el amarillo, rojo o azul, recuerden que se debe escribir en ingles; Después se definen las escalas de la lista de valores en los ejes X y Y de cada variable, en nuestro gráfico nosotros dejaremos las variables con una relación libre "free", puesto que no las vamos a utilizar, si quisiéramos que fueran puestas en el diagrama, tendríamos que darle " draw=TRUE", el problema es que como tenemos una gran cantidad de datos los nombres de cada uno de los datos del eje Y se superpondrán en la gráfica y no se entenderá nada. Finalmente lo que hacemos es modificar los puntos que utilizaremos dentro de las gráficas, con el argumento "col=1" le decimos a R el numero de color determinado que queremos utilizar en nuestro gráfico para los puntos del gráfico, nosotros utilizamos en color negro que es codificado con el numero 1, si lo quisiéramos de color rojo tendríamos que darle el numero 2 y si quisiéramos de color verde le daríamos 3 y asi sucesivamente, despues lo que hacemos sera definir el tamaño de los puntos que utilicemos en la grafica, nosotros utilizaremos un tamaño de 0.5 "cex= 0.5" pero eso depende del gusto de cada quien, tambien podemos cambiar el simbolo que queremos utilizar, por ejemplo pueden ponerse asteriscos, estrellas, triangulos o simplemente puntos como lo hicimos aqui "pch=16", si hubiésemos querido poner estrellas en vez de puntos debíamos darle "pch=11", se pueden encontrar todos los valores para cada simbolo que utilicemos en la pagina de quick R. Ya lo ultimo es poner las etiquetas de los dos ejes X y Y (xlab = "valor de la variable", ylab = "orden de los datos en el archivo texto") .
Ahora como había prometido, proporcionare un método no gráfico entre los cientos de miles que hay para encontrar outliers (este parece complicado pero no lo es) debido a que quiero matar dos pájaros de un solo tiro. En primer lugar entregare una función para la búsqueda de lo outliers, y en segundo lugar quiero dar la primer noción acerca de programación de funciones en R, mostrando como esta construida la función para la búsqueda de outliers explicando algunos argumentos de la función.
La función utilizada fue extraída y modificada de un ejercicio del libro "The art of R programming", exactamente el ejemplo 3.3.2 de la pagina 72, en el capitulo numero 3 (Matrices and arrays).
>buscar_outliers<- function (x){
buscar <-function (xcol){
media <-median(xcol)
desviacion<-abs(xcol-media)
return (which.max(desviacion))
}
return (apply(x, 2,buscar))
}
#(((!!!la función numero 1 es resaltada en celeste, y la función numero 2 en rojo que esta escrita dentro de la funcion 1, por ello tiene como fondo color celeste!!!)))Teniendo en cuenta que son dos funciones en una sola, lo primero que hacemos es darle un nombre a la función 1, por lo tanto llamaremos nuestra función 1 "buscar_outliers", esta función sera aplicada para la matriz "X" (buscar_outliers<- function (x)) que puede ser cualquier matriz, y lo que hace es aplicar la función 2 a todas las columnas de la matriz (x), para esto se utiliza el argumento "return (apply (x,2,buscar))" que indica lo que queremos obtener como resultado (return) de la función que en este caso son las posiciones de los outliers en cada una de las variables, por ello le indicamos que busque los outliers de las columnas/variables con el numero "2", si quisieramos buscar los outliers del las filas tendriamos que darle el numero "1", finalmente la palabra buscar es el nombre de la función 2 e indica que se aplique la funcion numero 2 a cada una de las columnas de la matriz x.
La función 2 esta escrita dentro de la función 1, lleva por nombre "buscar" como lo dije anteriormente, esta se aplica a las columnas de una matriz "X" (buscar<-function (xcol)), lo primero que hace es hallar la media de los valores en la variable ( media<-median(xcol)), después halla la desviación de cada uno de los valores absolutos de la columna restando a cada uno de los valores el valor de la media (desviacion <-abs(xcol-media)) y generando como resultado el valor máximo de desviación de las columnas (return(wich.max(desviacion))).
# De este modo, lo único que tenemos que dar es el nombre de nuestra función creada y especificar a cual matriz queremos aplicarla que en nuestro caso es "Sparrows", y nuestros resultados nos dirán cual es el dato en cada una de las variables que posee un valor outlier (por ejemplo para "wingcrd" el dato numero 536 debe ser considerado un outlier):
>buscar_outliers (Sparrows)
wingcrd flatwing tarsus head culmen nalospi
536 536 646 366 142 1090
wt bandstat initials Year Month Day
1 48 76 1 2 365
Location SpeciesCode Sex Age
6 2 1 1
Esto que acabo de presentar es un pequeñisimo ejemplo de como construir funciones, pero en la vida real no es necesario construir funciones para encontrar outliers (aunque es necesario construir muchísimas funciones para nuestras necesidades particulares), es importante empezar con algo bastante sencillo.
De hecho, existen bastantes paquetes disponibles en R para buscar outliers, PARA MI uno muy bueno y fácil de utilizar es el paquete "outliers", daré un par de ejemplos de las funciones disponibles es este paquete aplicándolas a nuestra matriz.
>install.packages ("outliers")
#Instalamos el paquete "outliers"
> library (outliers)
#Cargamos el paquete en nuestra área de trabajo.
>chisq.out.test(Sparrows$wingcrd)
# Realizo una búsqueda de outlier por medio del método de chi-cuadrado, este nos indicara el valor de la variable que debe ser tomada como un outlier en nuestra variable "wingcrd" de nuestra matriz "Sparrows"
(Para mas info acerca de este método...)
chi-squared test for outlier
data: Sparrows$wingcrd
X-squared = 19.7323, p-value = 8.908e-06
alternative hypothesis: highest value 68 is an outlier
>outlier(Sparrows)
#Busco los valores que deben ser considerados outliers en cada una de las variables de nuestra matriz "Sparrows". wingcrd flatwing tarsus head culmen nalospi
68.0 68.5 31.8 21.1 8.9 18.8
wt bandstat initials Year Month Day
9.5 2.0 9.0 2004.0 10.0 31.0
Location SpeciesCode Sex Age
1.0 3.0 0.0 2.0
2. Prueba de homogeneidad de varianzas.
Tal como los autores mencionan en el articulo y como es bien conocido, la homogeneidad de varianzas es una condición muy importante para realizar análisis de varianza (ANOVA), técnicas de regresión y análisis discriminante, entre otras. La homogeneidad de varianzas o homocedasticidad se presenta en un conjunto de datos en los cuales las variables endógenas presentan la misma varianza o muy cercana a la misma. (mas información).
El método gráfico que utilizaron los autores fue el método de "cajas y bigotes condicionales" (conditional boxplots) y este es el que voy a mostrar en el post ademas de mostrar una pequeñísima linea de comandos para realizar el test de Levene para homogeneidad de varianzas.
Para esta prueba utilizaremos la base de datos de los autores, esta contiene las tasas de consumo de alimento de una ave migradora (Limosa haemastica) y se quiere saber si el consumo de alimento varia entre sexos, periodos de tiempo o por una combinación de estas dos variables. Por lo tanto se debe asumir que la (i) variacion en las obseraciones delos sexos es similar, (ii) que la variación de las observaciones en diferentes periodos de tiempo es similar y (iii) que la variación entre periodos de tiempo y sexos es similar.
> library (lattice)
#Se carga el paquete "lattice" a nuestra área de trabajo.
#Llamamos a R la base de datos para hacer nuestro análisis, y le ponemos por nombre "limosa", el archivo en el que están guardado los datos fue llamado por los autores "Godwits.txt" y le damos "header=TRUE" para indicarle que la primera fila son los nombres de las variables.
>limosa$fSEX <- factor(limosa$SEX, levels = c(0, 1, 2),labels = c("No", "Hembras", "Machos"))
# Debido a que los datos están organizados de forma codificada, debemos crear un vector de nombres para los estados en la tabla; es decir, en el archivo para la variable "SEX" el estado 0 significa que tiene sexo desconocido, 1 que es una hembra y 2 que el individuo es un macho. Por ello se le dice a R que construya un factor "SEX" de la base de datos "limosa" (limosa$fSEX) a partir de los estados 0,1 y 2 de la variable "SEX" de la base "limosa" (limosa$SEX, levels = c(0, 1, 2)) y que etiquete cada uno de los estados respectivamente con los nombres "No", "Hembras" y "Machos" (labels = c("No", "Hembras", "Machos").
>limosa$fPERIOD <- factor(limosa$PERIOD, levels = c(0, 1, 2),labels = c("Verano", "Pre-migracion", "Invierno"))
#Se hace lo mismo que en el paso anterior, solo que esta vez lo que codificaremos es los tres periodos de tiempo con los nombres "Verano", "Pre-migracion" e "Invierno".
>grafico_homogeneidad<-bwplot(mgconsumed ~ fPERIOD | fSEX, data = limosa,
strip = strip.custom(bg = 'red'), subset = SEX!=0,
cex = 1.0, layout = c(2, 1),
xlab = "Periodo de migracion", ylab = "Tasa de consumo de alimento",
par.settings = list(
box.rectangle = list(col = 3),
box.umbrella = list(col = 2),
plot.symbol = list(cex = 1.5, col = 1)),
scales = list(x = list(relation = "same"),
y = list(relation = "same")))
# Lo primero que hacemos es darle un nombre a nuestro gráfico, seré "muy creativo" y le daré el nombre de "grafico_homogeneidad", después con una de las funciones de lattice (bwplot) creamos los gráficos de cajas y bigotes condicionales, para ello utilizamos la variable "mgconsumed" (mg de alimento consumido), le digo que tome los dos factores (periodo y sexo) y que uno depende del otro, con el argumento (fPERIOD | fSEX) le decimos a lattice que nos haga un gráfico por periodo de migración (es decir 3 gráficos) pero con la condición de que nos haga cada gráfico por periodo dependiendo del sexo (es decir 2 sexos X 3 periodos cada uno= 6 gráficos). Si le diéramos (fSEX | fPERIOD) nos haría el mismo numero de gráficos pero esta vez serian 3 periodos X 2 sexos cada uno y la configuración del gráfico seria distinta. Con las opciones del argumento "strip" lo que hacemos es modificar el cuadro de las gráficas y la organización de estas, para ello editamos con "strip.custom" la barra de etiquetas; primero le decimos que el color del fondo (background=bg), donde van las etiquetas de los sexos sea el color rojo (bg= 'red'), pero cabe aclarar que podriamos utilizar cualquier color que queramos, pero de debe poner en ingles el nombre del color!!. Después de esto, lo que hacemos es quitar el estado de sexo desconocido (0), para ello utilizamos el argumento (subset = SEX!=0) que indica a lattice que de los tres estados o subsets nos elimine el estado 0. Si deseamos incluir el estado de sexo "desconocido" (0) simplemente no utilizamos ese argumento.
Por otra parte le decimos a lattice el tamaño del punto que utilizaremos para indicar la media de la variable en cada caja con el argumento (cex = 1.0) que indica que el tamaño es de 1.0, pero nosotros podemos poner el tamaño que queramos como por ejemplo 1.5 o 0.4, seguido a esto lo que hacemos es decirle como organice cada uno de los sub-gráficos en el gráfico total, por ejemplo, podemos decirle que organice todos los sub-gráficos unos seguido de otros horizontalmente o unos encima de otros verticalmente, para eso lo que tenemos que decirle es en cuantas columnas y filas nos organice los sub-gráficos; En nuestro gráfico total le decimos que queremos una sola hilera o fila de gráficos con dos columnas, para ello utilizamos el comando layout = c(2, 1) (2 columnas y 1 fila), si quisiéramos que los sub-gráficos quedaran uno encima del otro le diríamos que los organice en una sola columna pero con dos filas (layout = c(1, 2)) y así sucesivamente dependiendo de cuantos sub-graficos tengamos.Finalmente lo que hacemos es la parte de edición del gráfico para que quede mas bonito o quede como nosotros queramos, para ello la linea de argumentos que utilizamos es la de ajuste de los parámetros del gráfico especificando con una lista los parámetros que queramos editar (par.settings =list(.....)). En nuestro gráfico solo editaremos 3 parámetros del gráfico (aunque uno puede editar muchísimos parámetros y solo haría falta adicionarlos a la lista), el primer parámetro es el color del rectángulo de las cajas, al cual le dimos un color verde que esta codificado con el numero "3", si quisiéramos aplicar el color rojo tendríamos que darle el numero "2" (box.rectangle = list(col = 3)) . El segundo parámetro que modificaremos es el color de el bigote de la caja (en ingles, la sombrilla de la caja), le damos el color numero "2" que es rojo (box.umbrella = list(col = 2)) . El tercer parámetro que vamos a editar es el del color de los puntos extremos (outliers) de las variables y el tamaño de estos en el gráfico, entonces le damos un tamaño de punto de 1.5 y elegimos el color negro para estos que esta codificado con el numero "1" (plot.symbol = list(cex = 1.5, col = 1)). Finalmente lo que hacemos es determinar la escala de los ejes pero en este caso el decimos que la relación de escalas entre los ejes es la misma (scales = list(x = list(relation = "same"), y = list(relation = "same"))) . Las etiquetas (labels) de los ejes fueron nombradas con el comando (xlab = "Periodo de migracion", ylab = "Tasa de consumo de alimento").
>grafico_homogeneidad
# Visualizamos nuestro gráfico.
Ahora bien, al igual que para encontrar outliers existen métodos no gráficos, para revisar sí tenemos homogeneidad de varianza en nuestros datos también hay muchos test disponibles en R, como la prueba de F que es la que utilizare para explicar su desarrollo en R, además del test de Levene entre otros.
> hembras<-subset(limosa, SEX==1, select = c(mgconsumed),header =T)
#Primero hacemos un substet que contenga la cantidad de alimento consumido por las hembras, utilizando la base de datos "limosa", el estado para la variable "SEX" utilizado para las hembras es el numero 1 (SEX==1), y seleccionamos los valores de cantidad de alimento consumido (select = c(mgconsumed)), correspondiente a los individuos que son hembras. Por ultimo le decimos que la primera fila de "limosa" es el nombre de las variables y también del nuevo subset (header =T).
> machos<-subset(limosa, SEX==2, select = c(mgconsumed),header =T)
#Hacemos lo propio, creando también el subset de los machos para comparar las varianzas.
>var.test(hembras$mgconsumed,machos$mgconsumed)
# Realizo mi test de homogeneidad de varianzas para los dos substes antes creados y obtenog mis resultados:
F test to compare two variances
data: hembras$mgconsumed and machos$mgconsumed
F = 0.6089, num df = 77, denom df = 127, p-value = 0.01891
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.4111176 0.9201650
sample estimates:
ratio of variances
0.60890273. Prueba de distribución normal de los datos.
Finalmente, para terminar este post tan largo, pero el cual creo que es muy útil (además es demasiado largo porque me esforcé muchísimo por explicar cada detalle minuciosamente, así que espero todo haya quedado muy claro), explicare la vía gráfica (histograma de frecuencias) para hacer la prueba de distribución normal en R, además de explicar una vía no gráfica entre una infinidad de formas de hacerlo en R.
>pericos <- read.table(file="SparrowsElphick.txt", header=TRUE)
#Para empezar lo primero que vamos a hacer es llamar nuestros datos (los datos de los autores, y crearemos un objeto en R llamado "pericos" (no son pericos, pero que hijuemadre, pongamole pericos), y le decimos que la primer fila de la tabla son los nombres de las variables con (header=T).
>pericos$fMonth<-factor(pericos$Month,
levels = c(5, 6, 7, 8, 9, 10),
labels = c("Mayo", "Junio", "Julio", "Agosto",
"Sept.", "Oct."))
#Despues creo un subset de datos, y para ello lo primero que debo hacer es escoger un factor que me separe algunas observaciones de otras, por lo tanto escojo la variable "mes" (month) del 5,6,7,8,9 y 10 que es como están codificados en la matriz y corresponden a los meses de mayo hasta Octubre) set de datos original, y extraigo las observaciones correspondientes de estos 6 meses.
>pericos$I1 <- pericos$fMonth =="Junio" |
pericos$fMonth =="Julio" |
pericos$fMonth =="Agosto" |
pericos$fMonth =="Mayo" |
pericos$fMonth =="Oct." |
pericos$fMonth =="Sept."
# Ahora creo un objeto en R con las observaciones obtenidas anteriormente, pero esta vez, cuando los mese ya han sido etiquetados y no están codificados numéricamente.
>hist(pericos$wt[pericos$I1],
xlab = "Peso (g)", breaks = 100,
main = "", ylab = "Frecuencia")
# Finalmente realizamos el histograma de frecuencias para visualizar si existe o no normalidad en los datos, utilizando los datos de peso de la matriz original (pericos$wt) (wt=weight=peso) de las observaciones correspondientes al factor antes creado que corresponde a algunos meses ([pericos$I1]), le ponemos nombre a los ejes con (xlab = "Peso (g) y "ylab = "Frecuencia"), y le decimos que el histogramas de frecuencias haga 100 rangos (breaks = 100), es decir 100 barritas de un peso en especifico cada una, pero podriamos decirle que haga el numero que nosotros deseemos.
De vez en cuando se deben dividir los datos, para así definir que parte esta cargando la mayor variabilidad o como se comportan los datos, por ejemplo en meses diferentes, para ello podemos hacer un histograma de frecuencias, para así revisar la normalidad de los datos de un mes en especifico.
>library(lattice)
#Cargo el paquete "lattice", para hacer graficos multiples.
>histogram( ~ wt | fMonth, type = "count",
xlab = "Peso (g)",
ylab = "Frecuencia",
nint=100,layout=c(3,2),
strip.left = strip.custom(bg = 'yellow'),
strip = F,
col.line = "red", col = "red",
scales = list(x = list(relation = "same"),
y = list(relation = "same"),
draw = TRUE),
subset = fMonth =="Mayo" |fMonth =="Junio" | fMonth == "Julio" |fMonth == "Agosto" |fMonth == "Sept." |fMonth == "Oct.",
data = pericos)
#Realizo mi histograma de frecuencias con los datos de peso (wt) de mi vector fMonth, y el tipo de histograma a realizar es del tipo "count", le pongo el nombre de los ejes de cada variable, con el argumento (nint=100), le digo que divida cada histograma en 100 rangos, con (layout=c(3,2), le digo que distribuya cada uno de los histogramas en 3 columnas y dos filas, los siguientes argumentos son para cambiar los colores de los gráficos las escalas de cada uno de los ejes para dejar las mismas para todos los gráficos; por ultimo especificamos el subset de datos que vamos a utilizar.
Y para realizar un aprueba de normalidad rápidamente en R, hay muchos caminos pero el que generalmente he usado yo es el test de Shapiro, disponible en algunos de mis posts anteriores.
uffffff
finalmente eso era todo lo de esta primera parte de exploración de datos para evitar errores estadísticos, y para los próximos dos posts les tendré los pasos restantes de todo el protocolo.
Proximo post:
Errores comunes en estadística: Exploración de datos (II) en R: Cantidad de ceros en los datos y Colinealidad entre las variables.
Suscribirse a:
Entradas (Atom)